© 2022 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

Published article:
S. Katayama and T. Ohtsuka, “Efficient Riccati recursion for optimal control problems with pure-state equality constraints,” 2022 American Control Conference (ACC), 2022, pp. 3579–3586, doi: 10.23919/ACC53348.2022.9867631.

Efficient Riccati recursion for optimal control problems
with pure-state equality constraints

Sotaro Katayama¹ and Toshiyuki Ohtsuka¹ ¹S. Katayama and T. Ohtsuka are with Department of System Science, Graduate School of Informatics, Kyoto University, Kyoto, Japan [email protected], [email protected]

Abstract

A novel approach to efficiently treat pure-state equality constraints in optimal control problems (OCPs) using a Riccati recursion algorithm is proposed. The proposed method transforms a pure-state equality constraint into a mixed state-control constraint such that the constraint is expressed by variables at a certain previous time stage. It is showed that if the solution satisfies the second-order sufficient conditions of the OCP with the transformed mixed state-control constraints, it is a local minimum of the OCP with the original pure-state constraints. A Riccati recursion algorithm is derived to solve the OCP using the transformed constraints with linear time complexity in the grid number of the horizon, in contrast to a previous approach that scales cubically with respect to the total dimension of the pure-state equality constraints. Numerical experiments on the whole-body optimal control of quadrupedal gaits that involve pure-state equality constraints owing to contact switches demonstrate the effectiveness of the proposed method over existing approaches.

I Introduction

Optimal control underlies the motion planning and control of dynamical systems such as trajectory optimization (TO) and model predictive control (MPC) [1]. TO achieves versatile and dynamically consistent planning by solving optimal control problems (OCPs). MPC leverages the same advantages as TO in real-time control by solving an OCP online within a particular sampling interval. It is essential, particularly for MPC, to solve direct OCPs within a short computational time, even if they involve highly complicated dynamics, a large dimensional state, and a long horizon.

Newton-type methods are the most practical methods used for solving OCPs in terms of the convergence speed. One of the most efficient algorithms that implement the Newton-type methods to solve both the single-shooting and multiple-shooting OCPs of large-scale systems is the Riccati recursion algorithm [1, 2]. The Riccati recursion algorithm scales only linearly with respect to the number of discretization grids of the horizon, in contrast to the direct methods (i.e., methods applying the Cholesky decomposition directly to the entire Hessian matrix) that scale cubically. For example, it was successfully applied to solve OCPs even for significantly complex systems such as legged robots with large degrees of freedom within very short computational times [3, 4].

However, there is a drawback of the Riccati recursion algorithm: it cannot efficiently treat pure-state equality constraints, which often arise, for example, in waypoint constraints, terminal constraints, switching constraints in hybrid systems such as legged robots [5, 6], and inequality constraints handled by active-set methods [7]. [7] extended the Riccati recursion algorithm for pure-state equality constraints and illustrated its effectiveness over the direct method for certain quadratic programming problems. However, the computational time of this method scales cubically with respect to the total dimension of pure-state equality constraints over the horizon, and it is inefficient when the total dimension is large. [8] proposed a projection approach to treat pure-state equality constraints with the Riccati recursion algorithm efficiently. However, this approach can only treat the equality constraints whose relative degree is 1, for example, velocity-level constraints of second-order systems, and cannot treat position-level constraints for such systems, which is a very common and practical class of constraints.

Other popular constraint-handling methods used with the Riccati recursion algorithm are the penalty function method and the augmented Lagrangian (AL) method. For example, [5] used the penalty function method and [6] used the AL method to treat pure-state equality constraints representing the switching constraints arising in OCPs involving quadrupedal gaits. A transformation of the linear-quadratic OCP into a dual problem for the efficient Riccati recursion [9] also used the penalty function method to treat the pure-state constraints. However, the penalty function method practically yields only the approximated solution, as illustrated by the numerical results obtained in [5]. The AL method can treat constraints better than the penalty function method, for example, it converges to the optimal solution even if the penalty parameter remains at a finite value. However, it generally lacks convergence speed compared with the Newton-type methods that achieve superlinear or quadratic convergence. The AL method essentially achieves linear convergence and it yields superlinear convergence only if the penalty parameter goes to infinity [10], which is an impractical assumption. For example, [6] used the AL method to consider the switching constraints in an OCP of quadruped bouncing motion. The AL method required a large number of the iterations (up to 300), although a simple 2D robot model was used and only one cycle of the bouncing motion was considered, which led to very low-dimensional (only four dimensions in all) pure-state equality constraints.

In this paper, we propose a novel approach to efficiently treat pure-state equality constraints in OCPs using a Riccati recursion algorithm. The proposed method transforms a pure-state equality constraint into a mixed state-control constraint such that the constraint is expressed by variables at a certain previous time stage. We show a relationship between an OCP with the original pure-state constraint (the original OCP) and an OCP with the transformed mixed state-control constraint (the transformed OCP); if the solution satisfies the first-order necessary conditions (FONC) and/or second-order sufficient conditions (SOSC) of the transformed OCP, then the solution also satisfies the FONC and/or SOSC of the original OCP. Therefore, if we find a solution that satisfies the SOSC of the transformed OCP, it is a local minimum of the original OCP. We then derive a Riccati recursion algorithm to solve the transformed OCP with linear time complexity in the grid number of the horizon, in contrast with the previous approach [7] that scales cubically with respect to the total dimension of pure-state equality constraints. Moreover, because the proposed method is in substance a Newton’s method for an optimization problem with equality constraints, the proposed method achieves superlinear or quadratic convergence, which distinguishes our approach from the penalty function method and the AL method in terms of convergence speed. We present numerical experiments of the whole-body optimal control of quadrupedal gaits that involve pure-state equality constraints owing to contact switches, which represent the position-level constraints of a second-order system, and demonstrate the effectiveness of the proposed method over the existing approaches, that is, the approach of [7] and the AL method.

This paper is organized as follows. In Section II, we transform an OCP with a pure-state equality constraint into an OCP with a mixed state-control equality constraint. In Section III, we derive a Riccati recursion algorithm to apply Newton’s method efficiently to the OCP with transformed mixed state-control equality constraints. In Section IV, the theoretical properties of the proposed transformation of OCPs are described. In Section V, the proposed method is compared with existing methods and its effectiveness is demonstrated in terms of computational time and convergence speed. In Section VI, we conclude with a brief summary and mention of future work.

Notation: We describe the partial derivatives of a differentiable function with respect to certain variables using a function with subscripts; that is, $f_{x}(x)$ denotes $\frac{\partial f}{\partial x}(x)$ and $g_{xy}(x,y)$ denotes $\frac{\partial^{2}g}{\partial x\partial y}(x,y)$ .

II Transformation of Optimal Control Problem with Pure-State Equality Constraints

II-A Original Optimal Control Problem

We consider the following discrete-time OCP: Find the state $x_{0},...,x_{N}\in\mathbb{R}^{n_{x}}$ and the control input $u_{0},...,u_{N-1}\in\mathbb{R}^{n_{u}}$ minimizing the cost function

J=\varphi(x_{N})+\sum_{i=0}^{N-1}L(x_{i},u_{i})

(1)

subject to the state equation

x_{i}+f(x_{i},u_{i})-x_{i+1}=0,\;\;i\in\left\{0,...,N-1\right\},

(2)

a pure-state equality constraint

\phi(x_{k})=0,\;\phi(x_{k})\in\mathbb{R}^{n_{c}},

(3)

and the initial state constraint

x_{0}-\bar{x}=0,\;\bar{x}\in\mathbb{R}^{n_{x}}.

(4)

In the following, we assume the form of the state as $x_{i}=\begin{bmatrix}q_{i}^{\rm T}&v_{i}^{\rm T}\end{bmatrix}^{\rm T}$ , where $q_{i}\in\mathbb{R}^{n}$ and $v_{i}\in\mathbb{R}^{n}$ represent the generalized coordinates and velocity of the system, respectively, and assume the form of the state equation as follows:

f(x_{i},u_{i}):=\begin{bmatrix}f^{(q)}(x_{i})\\ f^{(v)}(x_{i},u_{i})\end{bmatrix},\;\;f^{(q)}(x_{i}),\;f^{(v)}(x_{i},u_{i})\in\mathbb{R}^{n}.

(5)

We also assume that $k\geq 2$ , $n_{u}\leq n$ , $n_{c}\leq n$ , and that the constraint (3) depends only on the generalized coordinate, that is, its Jacobian is expressed as follows:

\phi_{x}(x_{k})=\begin{bmatrix}\phi_{q}(q_{k})&O\end{bmatrix},\;\;\phi_{q}(q_{k})\in\mathbb{R}^{n_{c}\times n}.

(6)

The state equation (5) mainly represents a second-order Lagrangian system with $n$ degrees of freedom, and a constraint (3), whose Jacobian is of form (6) represents a position-level constraint (that is, the relative degree of the constraint with respect to the control input is 2), which is a very common and practical class of problem settings.

II-B Transformation of Optimal Control Problem

To solve the aforementioned OCP efficiently, we transform the original pure-state equality constraint (3) into a mixed state-control equality constraint that is equivalent to (3) as long as (2) is satisfied. If (2) is satisfied,

	$\displaystyle x_{k}=\;$	$\displaystyle x_{k-2}+f(x_{k-2},u_{k-2})$
		$\displaystyle+f(x_{k-2}+f(x_{k-2},u_{k-2}),u_{k-1})$

holds and therefore

	$\displaystyle\phi($	$\displaystyle x_{k-2}+f(x_{k-2},u_{k-2})$
		$\displaystyle+f(x_{k-2}+f(x_{k-2},u_{k-2}),u_{k-1}))=0$		(7)

is equivalent to (3). Furthermore, because $\phi(\cdot)$ only depends on the generalized coordinate, (II-B) is equivalent to

\phi(x_{k-2}+f(x_{k-2},u_{k-2})+g(x_{k-2}+f(x_{k-2},u_{k-2})))=0,

(8)

where we define

g(x):=\begin{bmatrix}f^{(q)}(x)\\ 0\end{bmatrix}.

(9)

Therefore, the constraint (8) is equivalent to (3) if (2) is satisfied. Therefore, we consider (8) instead of (3) in the following. We herein summarize the original and transformed OCPs as follows:

Problem 1 – Original OCP: Find the solution $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ minimizing (1) subject to (2)–(4).

Problem 2 - Transformed OCP: Find the solution $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ minimizing (1) subject to (2), (4), and (8).

In fact, we have the following relations between the transformed OCP and the original OCP: If the solution $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ satisfies the FONC of the transformed OCP, it also satisfies the FONC of the original OCP. If the solution $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ satisfies the SOSC of the transformed OCP, it also satisfies the SOSC of the original OCP. Therefore, if we find the solution $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ that satisfies the SOSC of the transformed OCP, it is a local minimum of the original OCP. We show these theoretical points later in Section IV.

It should be noted that it is trivial to apply the proposed approach to constraints of relative degree 1: we transform (3) into a mixed state-control constraint represented by $x_{k-1}$ and $u_{k-1}$ in the same manner as mentioned before. Therefore, our approach comprises the approach of [8] that also involves constraint transformation but can treat only constraints of relative degree 1. The difference between our approach and that of [8] is that we introduce the constraint transformation in the original nonlinear problem and leverage the structure of the system (2), whereas in approach [8], constraint transformation is introduced in the linear subproblem arising in the Newton-type iterations without any assumptions in the state equation. As a result, only our approach can treat the practically important constraints of relative degree 2 with a theoretical justification for the transformation.

It should also be noted that the proposed approach is completely different from the classical transformation of pure-state equality constraints in continuous-time OCPs by considering their derivatives with respect to time, for example, in Section 3.4 of [11]. To explain this difference, we consider that there is a pure-state constraint of the form of (3) over a time interval. The classical method then transforms the pure-state constraint over the interval into a combination of the pure-state equality constraints (3) and $\frac{d}{dt}\phi(x)=0$ at a point on the interval and the mixed state-control constraint $\frac{d^{2}}{dt^{2}}\phi(x)=0$ over the interval, where $\frac{d}{dt}$ and $\frac{d^{2}}{dt^{2}}$ yield a kind of Lie derivative. Therefore, the classical method still needs to consider the pure-state equality constraints, whereas our approach transforms all pure-state equality constraints into the corresponding mixed state-control constraints.

II-C Optimality Conditions

We derive the optimality conditions, known as FONC, of the transformed OCP. We first define the Hamiltonian

H(x,u,\lambda):=L(x,u)+\lambda^{\rm T}f(x,u)

and

	$\displaystyle\tilde{H}(x,u,\lambda,\nu)$
	$\displaystyle:=H(x,u,\lambda)+\nu^{\rm T}\phi(x+f(x,u)+g(x+f(x,u))).$

We also define the intermediate time stages in which the constraint is not active as $\bar{I}:=\left\{1,...,k-3,k-1,...,N-1\right\}$ . The optimality conditions are then derived as follows [11]:

\varphi_{x}^{\rm T}(x_{N})-\lambda_{N}=0,

(10)

H_{x}^{\rm T}(x_{i},u_{i},\lambda_{i+1})+\lambda_{i+1}-\lambda_{i}=0

(11)

and

H_{u}^{\rm T}(x_{i},u_{i},\lambda_{i+1})=0

(12)

for $i\in\bar{I}$ ,

		$\displaystyle\tilde{H}_{x}^{\rm T}(x_{k-2},u_{k-2},\lambda_{k-1},\nu)+\lambda_{k-1}-\lambda_{k-2}$
		$\displaystyle={H}_{x}^{\rm T}(x_{k-2},u_{k-2},\lambda_{k-1})+\lambda_{k-1}-\lambda_{k-2}$
		$\displaystyle\;\;\;\;\;+(I+f_{x}^{\rm T}(x_{k-2},u_{k-2}))(I+g_{x}^{\rm T})\phi_{x}^{\rm T}\nu=0,$		(13)

and

		$\displaystyle\tilde{H}_{u}^{\rm T}(x_{k-2},u_{k-2},\lambda_{k-1},\nu)$
		$\displaystyle={H}_{u}^{\rm T}(x_{k-2},u_{k-2},\lambda_{k-1})$
		$\displaystyle\;\;\;\;\;+f_{u}^{\rm T}(x_{k-2},u_{k-2})(I+g_{x}^{\rm T})\phi_{x}^{\rm T}\nu=0,$		(14)

where $\lambda_{0},...,\lambda_{N}$ are the Lagrange multipliers with respect to (4) and (2), and $\nu$ is that with respect to (8). It should be noted that we have omitted the arguments from $\phi_{x}$ and $g_{x}$ in (II-C) and (II-C).

III Riccati Recursion

III-A Linearization for Newton’s Method

To apply Newton’s method for the transformed OCP, we linearize the optimality conditions. It should be noted that we adopt the direct multiple shooting method [12], that is, we consider $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ , $\lambda_{0},...,\lambda_{N}$ , and $\nu$ as the optimization variables.

III-A1 Terminal stage

At the terminal stage ( $i=N$ ), we have

Q_{xx,N}\Delta x_{N}-\Delta\lambda_{N}+\bar{l}_{x,N}=0,

(15)

where we define $Q_{xx,N}:=\varphi_{xx}(x_{N})$ . Further, we define $\bar{l}_{x,N}$ using the left-hand side of (10).

III-A2 Intermediate stages without equality constraint

In the intermediate stages without an equality constraint ( $i\in\bar{I}$ ), we have

A_{i}\Delta x_{i}+B_{i}\Delta u_{i}-\Delta x_{i+1}+\bar{x}_{i}=0,

(16)

Q_{xx,i}\Delta x_{i}+Q_{xu,i}\Delta u_{i}+A_{i}^{\rm T}\Delta\lambda_{i+1}-\Delta\lambda_{i}+\bar{l}_{x,i}=0,

(17)

and

Q_{xu,i}^{\rm T}\Delta x_{i}+Q_{uu,i}\Delta u_{i}+B_{i}^{\rm T}\Delta\lambda_{i+1}+\bar{l}_{u,i}=0,

(18)

where we define $A_{i}:=I+f_{x}(x_{i},u_{i})$ , $B_{i}\allowbreak:=f_{u}(x_{i},u_{i})$ , $Q_{xx,i}\allowbreak:=H_{xx}(x_{i},u_{i},\lambda_{i+1})$ , $Q_{xu,i}\allowbreak:=H_{xu}(x_{i},u_{i},\lambda_{i+1})$ , $Q_{uu,i}\allowbreak:=H_{uu}(x_{i},u_{i},\lambda_{i+1})$ . Further, we define $\bar{x}_{i}$ , $\bar{l}_{x,i}$ , and $\bar{l}_{u,i}$ using the left-hand sides of (2), (11), and (12), respectively.

III-A3 Intermediate stage with an equality constraint

At the intermediate stages with an equality constraint ( $i=k-2$ ), we have (16) for $i=k-2$ and have

	$\displaystyle Q_{xx,k-2}\Delta x_{k-2}$	$\displaystyle+Q_{xu,k-2}\Delta u_{k-2}+A_{k-2}^{\rm T}\Delta\lambda_{k-1}$
		$\displaystyle-\Delta\lambda_{k-2}+C^{\rm T}\Delta\nu+\bar{l}_{x,k-2}=0,$		(19)

	$\displaystyle Q_{xu,k-2}^{\rm T}\Delta x_{k-2}+Q_{uu,k-2}\Delta u_{k-2}+B_{k-2}^{\rm T}\Delta\lambda_{k-1}$
	$\displaystyle+D^{\rm T}\Delta\nu+\bar{l}_{u,k-2}=0,$		(20)

and

C\Delta x_{k-2}+D\Delta u_{k-2}+\bar{\phi}=0,

(21)

where we define $Q_{xx,k-2}\allowbreak:=\tilde{H}_{xx}(x_{k-2},u_{k-2},\lambda_{k-1},\nu)$ , $Q_{xu,k-2}\allowbreak:=\tilde{H}_{xu}(x_{k-2},u_{k-2},\lambda_{k-1},\nu)$ , $Q_{uu,k-2}\allowbreak:=\tilde{H}_{uu}(x_{k-2},u_{k-2},\lambda_{k-1},\nu)$ , $C\allowbreak:=\phi_{x}(I+g_{x})A_{k-2}$ , $D\allowbreak:=\phi_{x}(I+g_{x})B_{k-2}$ . We further define $\bar{l}_{x,k-2}$ , $\bar{l}_{u,k-2}$ , and $\bar{\phi}$ using the left-hand sides of (II-C), (II-C), and (8), respectively.

III-A4 Initial stage

Finally, we have

\Delta x_{0}+x_{0}-\bar{x}=0.

(22)

It should be noted that we can apply the Gauss-Newton Hessian approximation, which improves the computational speed when the constraints (2) and (8) are too complicated for their second-order derivatives to be computed. $Q_{xx,N}$ and $Q_{xx,i}$ , $Q_{xu,i}$ , and $Q_{uu,i}$ for $i\in\left\{0,...,N-1\right\}$ are then approximated using only the cost function (1) and do not depend on the Lagrange multipliers.

III-B Derivation of Riccati Recursion

We derive a Riccati recursion algorithm to solve the linear equations for Newton’s method (15)–(22) efficiently. As the standard Riccati recursion algorithm ([1, 2]), our goal is to derive a series of matrices $P_{i}$ and vectors $s_{i}$ such that

\Delta\lambda_{i}=P_{i}\Delta x_{i}-s_{i}

(23)

holds.

III-B1 Terminal stage

At the terminal stage ( $i=N$ ), we have

P_{N}=Q_{xx,N},\;s_{N}=-\bar{l}_{N}.

(24)

In the forward recursion, we have $\Delta x_{N}$ , and we compute $\Delta\lambda_{N}$ from (23).

III-B2 Intermediate stages without an equality constraint

At the intermediate stages without an equality constraint ( $i\in\bar{I}$ ), we have the following standard backward Riccati recursion ([1, 2]) for given $P_{i+1}$ and $s_{i+1}$ satisfying (23):

F_{i}:=Q_{xx,i}+A_{i}^{\rm T}P_{i+1}A_{i},

(25)

H_{i}:=Q_{xu,i}+A_{i}^{\rm T}P_{i+1}B_{i},

(26)

G_{i}:=Q_{uu,i}+B_{i}^{\rm T}P_{i+1}B_{i},

(27)

K_{i}:=-G_{i}^{-1}H_{i}^{\rm T},\;k_{i}:=-G_{i}^{-1}(B_{i}^{\rm T}P_{i+1}\bar{x}_{i}-B_{i}^{\rm T}s_{i+1}+\bar{l}_{u,i}),

(28)

and

P_{i}:=F_{i}-K_{i}^{\rm T}G_{i}K_{i},\;s_{i}:=A_{i}^{\rm T}(s_{i+1}-P_{i+1}\bar{x}_{i})-\bar{l}_{x,i}-H_{i}k_{i}.

(29)

In the forward recursion, for a particular value of $\Delta x_{i}$ , we compute $\Delta u_{i}$ from

\Delta u_{i}=K_{i}\Delta x_{i}+k_{i},

(30)

$\Delta\lambda_{i}$ from (23), and $\Delta x_{i+1}$ from (16).

III-B3 Intermediate stage with an equality constraint

At the intermediate stage with an equality constraint ( $i=k-2$ ), we first define (25)–(29) for $i=k-2$ for the specific values of $P_{k-1}$ and $s_{k-1}$ that satisfies (23). We then have the relations that are used in the forward recursion for $k-1$ as follows:

\begin{bmatrix}\Delta u_{k-2}\\ \Delta\nu\end{bmatrix}=\begin{bmatrix}K_{k-2}\\ M\end{bmatrix}\Delta x_{k-2}+\begin{bmatrix}k_{k-2}\\ m\end{bmatrix},

(31)

where we define

\begin{bmatrix}K_{k-2}\\ M\end{bmatrix}:=-\begin{bmatrix}G_{k-2}&D^{\rm T}\\ D&O\end{bmatrix}^{-1}\begin{bmatrix}H_{k-2}^{\rm T}\\ C\end{bmatrix}

(32)

and

	$\displaystyle\begin{bmatrix}k_{k-2}\\ m\end{bmatrix}:=-$	$\displaystyle\begin{bmatrix}G_{k-2}&D^{\rm T}\\ D&O\end{bmatrix}^{-1}$
		$\displaystyle\begin{bmatrix}B_{k-2}^{\rm T}P_{k-1}\bar{x}_{k-2}-B_{k-2}^{\rm T}s_{k-1}+\bar{l}_{u,k-2}\\ \bar{\phi}\end{bmatrix}.$		(33)

We then obtain the backward recursions

P_{k-2}:=F_{k-2}-\begin{bmatrix}K_{k-2}^{\rm T}&M^{\rm T}\end{bmatrix}\begin{bmatrix}G_{k-2}&D^{\rm T}\\ D&O\end{bmatrix}\begin{bmatrix}K_{k-2}\\ M\end{bmatrix}

(34)

and

	$\displaystyle s_{k-2}:=\;$	$\displaystyle A_{k-2}^{\rm T}(s_{k-1}-P_{k-1}\bar{x}_{k-2})-\bar{l}_{x,k-2}-H_{k-2}k_{k-2}$
		$\displaystyle-C^{\rm T}m.$		(35)

III-C Algorithm, Convergence, and Computational Analysis

We summarize the single Newton iteration, that is, the computation of the Newton direction for a particular solution, using the proposed Riccati recursion algorithm (Algorithm 1). In the first step, we formulate the linear equations of Newton’s method, that is, we compute the coefficient matrices and residuals of (15)–(22) (line 1). This step can leverage parallel computing. Second, we perform the backward Riccati recursion and compute $P_{i}$ and $z_{i}$ for $i\in\left\{0,...,N\right\}$ (lines 4–11). Third, we perform the forward Riccati recursion and compute the Newton directions for all the variables (lines 12–20).

Because the proposed method is in substance a Newton’s method for an optimization problem with equality constraints, it achieves superlinear or quadratic convergence, for example, by Proposition 4.4.3 of [10], which distinguishes the convergence behavior of the proposed method from that of the AL method, popularly used to treat the pure-state equality constraints with the Riccati recursion algorithm. The AL method achieves superlinear convergence only if the penalty parameter goes to infinity, which is an impractical assumption; otherwise, its convergence rate is just linear.

It should be noted that we can trivially apply the proposed method to OCPs with multiple pure-state equality constraints on the horizon. When there are multiple time stages involving constraint (3) on the horizon, we compute the coefficient matrices and residuals in (III-A3)–(21) for each of the time stages in line 1 of Algorithm 1, apply line 7 of Algorithm 1 for each of the time stages in the backward Riccati recursion, and apply line 15 of Algorithm 1 for each of the time stages in the forward Riccati recursion.

The proposed method is particularly efficient when there are several stages with pure-state equality constraints on the horizon. Suppose that there is an $n_{c,i}$ -dimensional pure-state equality constraint at each time stage $i$ of the horizon ( $n_{c,i}=0$ if there is no constraint at stage $i$ ). The proposed method then computes the inverse of a matrix whose size is $(n_{u}+n_{c,i})\times(n_{u}+n_{c,i})$ at each time stage in the backward recursion. In contrast, the previous approach of [7] requires the computation of the inverse of a matrix of size $(\sum_{i=0}^{N}n_{c,i})\times(\sum_{i=0}^{N}n_{c,i})$ . Broadly speaking, the computational burden of the proposed method with respect to the grid number $N$ is $O(N)$ , whereas that of the approach in [7] is $O(N^{3})$ .

It should be noted that it is easy to apply the proposed method to the single-shooting methods, which are popular in robotic applications [3, 4, 5, 6], by considering only $u_{0},...,u_{N-1}$ and $\nu$ as the decision variables. In the single-shooting methods, before line 1 of Algorithm 1, we first compute $x_{0},...,x_{N}$ based on $x(t_{0})$ and $u_{0},...,u_{N-1}$ using the state equation (2) sequentially. Further, we compute $\lambda_{N},...,\lambda_{0}$ using (10), (11), and (II-C) based on $u_{0},...,u_{N-1}$ , $\nu$ , and $x_{0},...,x_{N}$ , respectively, in the backward recursion (lines 5–11 of Algorithm 1). We can then compute the Newton directions $\Delta u_{0},...,\Delta u_{N-1}$ and $\Delta\nu$ using the same (or a similar) forward recursion (lines 12–20 of Algorithm 1).

Algorithm 1 Computation of Newton direction using proposed Riccati recursion

1:Initial state

{x}(t_{0})

, the current solution

x_{0},...,x_{N}

u_{0},...,u_{N-1}

, and Lagrange multipliers

\lambda_{0},...,\lambda_{N}

\nu

2:Newton directions

\Delta x_{0}

, …,

\Delta x_{N}

\Delta u_{0}

, …,

\Delta u_{N-1}

\Delta\lambda_{0}

, …,

\Delta\lambda_{N}

, and

\Delta\nu

3:for

i=0,\cdots,N

do in parallel

4: Computes the matrices and residuals in (15)–(22).

5:end for

6:Compute

P_{N}

and

z_{N}

from (24).

7:for

i=N,\cdots,0

do in serial

8: if

i=k-2

then

9: Computes

P_{k-2}

and

z_{k-2}

from (25)–(27), (32), (III-B3), and (III-B3).

10: else

11: Computes

P_{i}

and

z_{i}

from (25)–(29).

12: end if

13:end for

14:Compute

\Delta x_{0}

from (22).

15:for

i=0,\cdots,N-1

do in serial

16: if

i=k-2

then

17: Compute

\Delta u_{k-2}

\Delta\nu

\Delta\lambda_{k-2}

, and

\Delta x_{k-1}

from (31), (23), and (16).

18: else

19: Compute

\Delta u_{i}

\Delta\lambda_{i}

, and

\Delta x_{i+1}

from (30), (23), and (16).

20: end if

21:end for

22:Compute

\Delta\lambda_{N}

from (23).

IV Theoretical Properties of Optimal Control Problem Transformation

We show the theoretical relationships between the transformed OCP and the original OCP in this section. The first theorem concerns a stationary point of the transformed OCP and a stationary point of the original OCP.

Theorem IV.1

Suppose that $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ , $\lambda_{0},...,\lambda_{N}$ , and $\nu$ satisfy the FONC of the transformed OCP. Then, there exist the Lagrange multipliers $\lambda_{0}^{*},...,\lambda_{N}^{*}$ and $\nu^{*}$ such that $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ , $\lambda_{0}^{*},...,\lambda_{N}^{*}$ , and $\nu^{*}$ satisfy the FONC of the original OCP.

Proof:

First, we define the intermediate time stages without the active constraints of the original OCP as $\tilde{I}:=\left\{0,...,k-1,k+1,...,N-1\right\}$ . The FONC of the original OCP is then expressed by (2)–(4), (10), and (11) for $i\in\tilde{I}$ ,

H_{x}^{\rm T}(x_{k},u_{k},\lambda_{k+1}^{*})+\phi_{x}^{\rm T}\nu^{*}+\lambda_{k+1}^{*}-\lambda_{k}^{*}=0,

(36)

and (12) for $i\in\left\{0,...,N-1\right\}$ , in which (2)–(4) are trivially satisfied because $x_{0},...,x_{N}$ and $u_{0},...,u_{N-1}$ satisfy the FONC of the transformed OCP. It should be noted that because $x_{0},...,x_{N}$ and $u_{0},...,u_{N-1}$ satisfy (2) and $\phi_{x}(\cdot)$ only depends on the generalized coordinate,

	$\displaystyle\phi_{x}(x_{k})=\phi_{x}($	$\displaystyle x_{k-2}+f(x_{k-2},u_{k-2})$
		$\displaystyle+g(x_{k-2}+f(x_{k-2},u_{k-2})))$		(37)

holds. Therefore, we describe both the left- and right-hand sides of (IV) as $\phi_{x}$ in this proof. Let

\lambda_{i}^{*}=\lambda_{i},\;\;i\in\left\{0,...,k-2,k+1,...,N\right\}.

(38)

Then, (10), (11) for $i\in\left\{0,...,k-3,k+1,...,N\right\}$ , and (12) for $i\in\left\{0,...,k-2,k+1,...,N\right\}$ of the original OCP are reduced to those of the transformed OCP and are, therefore, satisfied. Furthermore, let

\nu^{*}=\nu,\;\lambda_{k}^{*}=\lambda_{k}+\phi_{x}^{\rm T}\nu^{*},\;\lambda_{k-1}^{*}=\lambda_{k-1}+(I+g_{x}^{\rm T})\phi_{x}^{\rm T}\nu^{*}.

(39)

Then, (11) and (12) for $i=k-2,k-1$ and (36) of the original OCP are also reduced to (II-C), (II-C), (11) and (12) for $i=k-1$ , and (11) for $i=k$ of the transformed OCP, respectively, noting that $\phi_{x}f_{x}(x_{k-1},u_{k-1})=\phi_{x}g_{x}$ and $\phi_{x}f_{u}(x_{k-1},u_{k-1})=O$ , and are, therefore, satisfied, which completes the proof. ∎

From the proof of Theorem 4.1, we can obtain the Lagrange multipliers at a stationary point of the original OCP corresponding to those of the transformed OCP. The following theorem concerns the sufficiency of the optimality.

Theorem IV.2

Suppose that $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ , $\lambda_{0},...,\lambda_{N}$ , and $\nu$ satisfy the SOSC of the transformed OCP. Then, there exist the Lagrange multipliers $\lambda_{0}^{*},...,\lambda_{N}^{*}$ , and $\nu^{*}$ such that $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ , $\lambda_{0}^{*},...,\lambda_{N}^{*}$ , and $\nu^{*}$ satisfy the SOSC of the original OCP.

Proof:

Because $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ , $\lambda_{0},...,\lambda_{N}$ , and $\nu$ also satisfy the FONC of the transformed OCP, we have $\lambda_{0}^{*},...,\lambda_{N}^{*}$ and $\nu^{*}$ defined by (38) and (39) such that $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ , $\lambda_{0}^{*},...,\lambda_{N}^{*}$ , and $\nu^{*}$ satisfy the FONC of the original OCP from Theorem 4.1. From the assumption of the SOSC of the transformed OCP, we have

\delta x_{N}^{\rm T}Q_{xx,N}\delta x_{N}+\sum_{i=0}^{N-1}\begin{bmatrix}\delta x_{i}^{\rm T}\\ \delta u_{i}^{\rm T}\end{bmatrix}^{\rm T}\begin{bmatrix}Q_{xx,i}&Q_{xu,i}\\ Q_{ux,i}^{\rm T}&Q_{uu,i}\end{bmatrix}\begin{bmatrix}\delta x_{i}\\ \delta u_{i}\end{bmatrix}>0

(40)

for arbitrary $\delta x_{i}$ and $\delta u_{i}$ satisfying $\delta x_{0}=0$ , $(I+f_{x,i})\delta x_{i}+f_{u,i}\delta u_{i}-\delta x_{i+1}=0$ for $i\in\left\{0,...,N-1\right\}$ and $\phi_{x}(I+g_{x})(I+f_{x,k-2})\delta x_{k-2}+\phi_{x}(I+g_{x})f_{u,k-2}\delta u_{k-2}=0$ , where we describe $f_{x,i}:=f_{x}(x_{i},u_{i})$ and $f_{u,i}:=f_{u}(x_{i},u_{i})$ for $i\in\left\{0,...,N-1\right\}$ . We introduce the Hessians of the original OCP as $Q_{xx,N}^{*}:=\varphi_{xx}(x_{N})=Q_{xx,N}$ and $Q_{xx,i}^{*}\allowbreak:=H_{xx}(x_{i},u_{i},\lambda_{i+1}^{*})$ , $Q_{xu,i}^{*}\allowbreak:=H_{xu}(x_{i},u_{i},\lambda_{i+1}^{*})$ , and $Q_{uu,i}^{*}\allowbreak:=H_{uu}(x_{i},u_{i},\lambda_{i+1}^{*})$ for $i\in\tilde{I}$ . Further, we introduce $Q_{xx,k}^{*}\allowbreak:=H_{xx}(x_{k},u_{k},\lambda_{k+1}^{*},\nu^{*})$ , $Q_{xu,k}^{*}\allowbreak:=H_{xu}(x_{k},u_{k},\lambda_{k+1}^{*},\nu^{*})$ , and $Q_{uu,k}^{*}\allowbreak:=H_{uu}(x_{k},u_{k},\lambda_{k+1}^{*},\nu^{*})$ . We can then complete the proof if

		$\displaystyle{\delta x_{N}^{}}^{\rm T}Q_{xx,N}^{}\delta x_{N}^{*}$
		$\displaystyle+\sum_{i=0}^{N-1}\begin{bmatrix}{\delta x_{i}^{}}^{\rm T}\\ {\delta u_{i}^{}}^{\rm T}\end{bmatrix}^{\rm T}\begin{bmatrix}Q_{xx,i}^{}&Q_{xu,i}^{}\\ {Q_{ux,i}^{}}^{\rm T}&Q_{uu,i}^{}\end{bmatrix}\begin{bmatrix}\delta x_{i}^{}\\ \delta u_{i}^{}\end{bmatrix}>0$		(41)

holds for arbitrary $\delta x_{i}^{*}$ and $\delta u_{i}^{*}$ satisfying $\delta x_{0}^{*}=0$ , $(I+f_{x,i})\delta x_{i}^{*}+f_{u,i}\delta u_{i}^{*}-\delta x_{i+1}^{*}=0$ for $i\in\left\{0,...,N-1\right\}$ and $\phi_{x}\delta x_{k}^{*}=0$ . First, we can see that the subspace of the feasible variations $\delta x_{i}^{*}$ and $\delta u_{i}^{*}$ is identical to that of $\delta x_{i}$ and $\delta u_{i}$ because (IV) holds. Then, we consider $\delta x_{i}^{*}$ and $\delta u_{i}^{*}$ as being identical to $\delta x_{i}$ and $\delta u_{i}$ . Next, by substituting (38) and (39) into the Hessians of the original OCP, we obtain $Q_{xx,i}^{*}=Q_{xx,i}$ , $Q_{xu,i}^{*}=Q_{xu,i}$ , and $Q_{uu,i}^{*}=Q_{uu,i}$ for $i\in\left\{0,...,k-3,k+1,...,N-1\right\}$ , $Q_{xx,k}^{*}=Q_{xx,k}+\nu\cdot\phi_{xx}$ , $Q_{xu,k}^{*}=Q_{xu,k}$ , $Q_{uu,k}^{*}=Q_{uu,k}$ , $Q_{xx,k-1}^{*}=Q_{xx,k-1}+(\phi_{x}^{\rm T}\nu)\cdot f_{xx,k-1}$ , $Q_{xu,k-1}^{*}=Q_{xu,k-1}+(\phi_{x}^{\rm T}\nu)\cdot f_{xu,k-1}$ , $Q_{uu,k-1}^{*}=Q_{uu,k-1}+(\phi_{x}^{\rm T}\nu)\cdot f_{uu,k-1}$ ,

	$\displaystyle Q_{xx,k-2}^{*}=Q_{xx,k-2}$
	$\displaystyle-(I+f_{x,k-2}^{\rm T})((\phi_{x}^{\rm T}\nu)\cdot g_{xx})(I+f_{x,k-2})$
	$\displaystyle-(I+f_{x,k-2}^{\rm T})(I+g_{x}^{\rm T})(\nu\cdot\phi_{xx})(I+g_{x})(I+f_{x,k-2}),$

	$\displaystyle Q_{xu,k-2}^{*}=Q_{xu,k-2}-(I+f_{x,k-2}^{\rm T})((\phi_{x}^{\rm T}\nu)\cdot g_{xx})f_{u,k-2}$
	$\displaystyle-(I+f_{x,k-2}^{\rm T})(I+g_{x}^{\rm T})(\nu\cdot\phi_{xx})(I+g_{x})f_{u,k-2},$

and

Q_{uu,k-2}^{*}=Q_{uu,k-2}-f_{u,k-2}^{\rm T}(I+g_{x}^{\rm T})(\nu\cdot\phi_{xx})(I+g_{x})f_{u,k-2},

where $f_{xx,i}:=f_{xx}(x_{i},u_{i})$ , $f_{xu,i}:=f_{xu}(x_{i},u_{i})$ , $f_{uu,i}:=f_{uu}(x_{i},u_{i})$ , and the notation “ $\cdot$ " denotes vector–tensor multiplication. Since the FONC of the original OCP holds and $\nu\cdot\phi_{xx}=\begin{bmatrix}\nu\cdot\phi_{qq}&O\\ O&O\end{bmatrix}$ , we have $(\nu\cdot\phi_{xx})(I+f_{x,k-1})=(\nu\cdot\phi_{xx})(I+g_{x})$ and $(\nu\cdot\phi_{xx})f_{u,k-1}=O$ , which yields $\delta x_{k}^{\rm T}(\nu\cdot\phi_{xx})\delta x_{k}=\delta x_{k-1}^{\rm T}(I+g_{x}^{\rm T})(\nu\cdot\phi_{xx})(I+g_{x})\delta x_{k-1}$ . In addition, from the structures of (5), (6), and (9), we have $(\phi_{x}^{\rm T}\nu)\cdot f_{xx,k-1}=(\phi_{x}^{\rm T}\nu)\cdot g_{xx}$ , $(\phi_{x}^{\rm T}\nu)\cdot f_{xu,k-1}=O$ , and $(\phi_{x}^{\rm T}\nu)\cdot f_{uu,k-1}=O$ . By substituting these relations, the left-hand side of (IV) is reduced to the left-hand side of (40), which completes the proof. ∎

We summarize the property of the proposed transformation in the next proposition:

Proposition IV.3

Suppose that $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ , $\lambda_{0},...,\lambda_{N}$ , and $\nu$ satisfy the SOSC of the transformed OCP. Then, the solution $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ is a strict local minimum of the original OCP.

Proof:

Because the solution $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ with the Lagrange multipliers satisfies the SOSC of the original OCP, as indicated by Theorem 4.2, the solution $x_{0},...,x_{N}$ , $u_{0},...,u_{N-1}$ is a strict local minimum of the original OCP. ∎

V Numerical Experiments on Whole-Body Quadrupedal Gaits Optimization

V-A Experimental Settings

To demonstrate the effectiveness of the proposed method over existing methods, we conducted numerical experiments on the whole-body optimal control of a quadrupedal robot ANYmal for various gaits. The equation of motion of the full 3D model of the quadrupedal robot is of the form of (5). Moreover, a pure-state constraint whose Jacobian is of the form (6) is imposed just before each impact between the leg and the ground, which is termed the switching constraint [5, 6]. We compare the following three Riccati recursion algorithms based on the direct-multiple shooting method and Gauss-Newton Hessian approximation with various constraint handling methods:

•

The proposed method
•

The Riccati recursion with pure-state constraints [7]
•

The AL method [13]

We implemented these three algorithms in C++ and used Pinocchio [14], an efficient C++ library used for rigid-body dynamics and its analytical derivatives, to compute the dynamics and its derivatives of the quadrupedal robot. We used OpenMP [15] for parallel computing (e.g., line 1 of Algorithm 1) and four threads through the following experiments. To consider the practical situation, we also imposed inequality constraints on the joint angle limits, joint angular velocity limits, and joint torque limits of each joint. We used the primal-dual interior point method [16] with fixed barrier parameters for the inequality constraints. None of the three methods used line search; they only used the fraction-to-boundary rule [16] for step-size selection. We fixed the instants of the impact between the robot and the ground in the following experiments and did not treat them as optimization variables as in [5, 6], to focus on the evaluation of the constraint-handling methods. All experiments were conducted on a laptop with a hexa-core CPU Intel Core i9-8950HK @2.90 GHz.

In the following two experiments, we considered that the OCP converges when the $l_{2}$ -norm of the residuals in the Karush—Kuhn–Tucker (KKT) conditions, which we refer to as the KKT error, becomes smaller than a prespecified threshold. The KKT conditions are composed of the FONC and primal and dual residuals in the inequality constraints. For example, the KKT conditions of the proposed method are composed of (2), (4), (8), (10)–(II-C), and the residuals in the inequality constraints. The KKT conditions of the Riccati recursion with pure-state constraints [7] and the AL method are only slightly different depending on the method used to treat pure-state equality constraints.

V-B Trotting Gait for Different Numbers of Steps

First, we evaluated the performances of the three methods for different total dimensions of pure-state equality constraints by considering the trotting gaits of ANYmal with different numbers of trotting steps. A six-dimensional (three-dimensional for each impact leg) pure-state equality constraint (switching constraint) was imposed on the OCPs for each trotting step. We chose the number of trotting steps from 2, 4, 6, 8, and 10 and measured the CPU time per Newton iteration and the total CPU time until convergence (we chose $1.0\times 10^{-10}$ as the convergence tolerance of the KKT errors). The settings used for the OCPs (horizon length $T$ , number of grids $N$ , and total dimension of equality constraints (3)) are listed in Table I. We carefully tuned the parameters of the AL method [13]. For example, we chose the initial penalty parameter as $p=5$ and the penalty update value as $\beta=8$ ; that is, the AL method updates the penalty parameter as $p\leftarrow\beta p$ when the KKT error excluding the constraint violation (3) is smaller than a tolerance that is also tuned carefully [13].

Figure 1 depicts the CPU time per Newton iteration (left figure) and the total CPU time until convergence (right figure) of each method. As depicted in the left figure of Fig. 1, the CPU time per Newton iteration in the proposed method was almost the same as that in the AL method, whereas the Riccati recursion with pure-state constraints [7] took more computational time when compared with the other two methods. The right figure of Fig. 1 also indicates that the proposed method achieved the fastest convergence. The proposed method took exactly the same number of iterations (approximately 20) until convergence as the Riccati recursion with pure-state constraints [7] in all the cases. Therefore, the proposed method was faster than it in terms of the total CPU time as in the case of the per Newton iteration. The AL method was significantly slower than the other two methods with respect to the total CPU time because it required approximately 80 iterations in all the cases, although we carefully tuned the AL algorithm.

TABLE I: Settings of OCPs for each number of trotting steps

No. of trotting steps	2	4	6	8	10
Horizon length $T$	1.55	2.55	3.55	4.55	5.55
No. of grids $N$	35	59	83	107	131
Total dim. of (3)	12	24	36	48	60

Refer to caption — Figure 1: (Left) CPU time per Newton iteration and (right) total CPU time until convergence for different numbers of trotting steps in the proposed method (Proposed), Riccati recursion with pure-state constraints [7] (pure state), and the AL method (AL).

V-C Trotting, Jumping, and Running Gait Problems

Next, we investigated the performances of the proposed method in three different quadrupedal gaits: trotting, jumping, and running gaits, of which the jumping and running gaits are particularly highly nonlinear and complicated problems. Each jumping step imposes a 12-dimensional pure-state equality constraint, and each running step imposes a 6-dimensional one. We summarize the settings of each problem (horizon length $T$ , number of grids $N$ , number of steps, total dimension of equality constraints (3), tolerance of convergence, and initial penalty parameter of the AL method $p_{\rm init}$ ) in Table II. As done in the preceding example, we carefully tuned the parameters of the AL method, that is, $p_{\rm init}$ and the update rule of the penalty and the Lagrange multiplier, for each problem. We measured the KKT errors with respect to the number of iterations, and the total number of iterations and CPU time until convergence.

Figure 2 depicts the $\log_{10}$ -scaled KKT error of each method for the three gait problems with respect to the number of iterations. We can see that the convergence behavior of the proposed method was almost the same as that of the Riccati recursion with pure-state constraints [7]. In contrast, the AL method resulted in significantly slow convergence because it needs to update the penalty parameter and Lagrange multiplier to reduce the constraint violation, which we can see in the peaks in the KKT error of the AL method in Fig. 2. Figure 3 indicates the number of iterations and total CPU time until convergence of each method. We can see that the total number of iterations of the proposed method was almost the same as that of the Riccati recursion with pure-state constraints [7], whereas the AL method required a significantly large number of iterations. In addition, as each iteration of the proposed method was faster than that of the Riccati recursion with pure-state constraints [7], as in the previous experiment, the proposed method achieved the fastest convergence.

TABLE II: Settings of OCPs for trotting, jumping, and running gaits.

Gait type	Trotting	Jumping	Running
Horizon length $T$	6.05	5	7
No. of grids $N$	143	107	346
No. of steps	11	3	26
Total dim. of (3)	66	36	156
KKT tolerance	$1.0\times 10^{-10}$	$1.0\times 10^{-10}$	$1.0\times 10^{-8}$
$p_{\rm init}$	5	1000	5

VI Conclusions

We proposed a novel approach to efficiently treat pure-state equality constraints in OCPs with a Riccati recursion algorithm. The proposed method transforms a pure-state equality constraint into a mixed state-control constraint such that the constraint is expressed by variables at a certain previous time stage. We derived a Riccati recursion algorithm to solve the transformed OCP with linear time complexity in the grid number of the horizon, in contrast to the previous approach [7], which scaled cubically with respect to the total dimension of the pure-state equality constraints. Because the proposed method is essentially a Newton’s method for an optimization problem with equality constraints, the proposed method achieves superlinear or quadratic convergence, which distinguishes our approach from the penalty function method and the AL method in terms of the convergence property. We showed that if the solution satisfies the FONC and/or SOSC of the transformed OCP, then the solution also satisfies the FONC and/or SOSC of the original OCP. Therefore, if we find a solution that satisfies the SOSC of the transformed OCP, it is a local minimum of the original OCP. We performed numerical experiments on the whole-body optimal control of quadrupedal gaits that involve pure-state equality constraints owing to contact switches and demonstrated the effectiveness of the proposed method over the approach of [7] and the AL method.

Our future work will include applying the proposed method with switching time optimization problems [5, 6]. We further extend the proposed method to constraints whose relative degree is larger than 2.

Acknowledgment

This work was partly supported by JST SPRING, Grant Number JPMJSP2110.

References

[1] K. Rawlings, D. Q. Mayne, and M. Diehl, Model Predictive Control: Theory, Computation, and Design. Nob Hill Publishing, LCC, 2017.
[2] G. Frison, “Algorithms and methods for high-performance model predictive control,” Ph.D. dissertation, Technical University of Denmark, 2016.
[3] J. Koenemann, A. Del Prete, Y. Tassa, E. Todorov, O. Stasse, M. Bennewitz, and N. Mansard, “Whole-body model-predictive control applied to the HRP-2 humanoid,” in 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2015, pp. 3346–3351.
[4] M. Neunert, M. Stäuble, M. Giftthaler, C. D. Bellicoso, J. Carius, C. Gehring, M. Hutter, and J. Buchli, “Whole-body nonlinear model predictive control through contacts for quadrupeds,” IEEE Robotics and Automation Letters, vol. 3, no. 3, pp. 1458–1465, 2018.
[5] F. Farshidian, M. Neunert, A. W. Winkler, G. Rey, and J. Buchli, “An efficient optimal planning and control framework for quadrupedal locomotion,” in 2017 IEEE International Conference on Robotics and Automation (ICRA), 2017, pp. 93–100.
[6] H. Li and P. M. Wensing, “Hybrid systems differential dynamic programming for whole-body motion planning of legged robots,” IEEE Robotics and Automation Letters, vol. 5, no. 4, pp. 5448–5455, 2020.
[7] A. Sideris and L. A. Rodriguez, “A Riccati approach for constrained linear quadratic optimal control,” International Journal of Control, vol. 84, no. 2, pp. 370–380, 2011.
[8] M. Giftthaler and J. Buchli, “A projection approach to equality constrained iterative linear quadratic optimal control,” in 2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids), 2017, pp. 61–66.
[9] D. Axehill, “Applications of integer quadratic programming in control and communication,” Ph.D. dissertation, Institutionen för systemteknik, 2005.
[10] D. P. Bertsekas, Nonlinear Programming, 3rd ed. Athena Scientific, 2016.
[11] A. E. Bryson and Y.-C. Ho, Applied Optimal Control: Optimization, Estimation, and Control. CRC Press, 1975.
[12] H. Bock and K. Plitt, “A multiple shooting algorithm for direct solution of optimal control problems,” in 9th IFAC World Congress, 1984, pp. 1603–1608.
[13] J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. Springer, 2006.
[14] J. Carpentier, G. Saurel, G. Buondonno, J. Mirabel, F. Lamiraux, O. Stasse, and N. Mansard, “The Pinocchio C++ library – A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives,” in International Symposium on System Integration (SII), 2019, pp. 614 – 619.
[15] L. Dagum and R. Menon, “OpenMP: An industry-standard API for shared-memory programming,” IEEE Computational Science & Engineering, vol. 5, no. 1, p. 46–55, 1998.
[16] A. Wächter and L. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematical Programming, vol. 106, pp. 25–57, 2006.