Optimal Control, Numerical Optimization, Interior Point Method

1 Introduction↩︎

We start by introducing the core concepts motivating our work.

1.1 Optimal Control↩︎

A discrete-time optimal control problem is an optimization problem of the form

\[\begin{align} &\min_{x_0, u_0, \ldots, x_N} &\sum \limits_{i=0}^{N-1} f_i(x_i, u_i) + f_N(x_N) \\ &s.t. & x_0 = s_0, \\ & & \forall i \in \lbrace 0, \ldots, N-1 \rbrace, x_{i+1} = d_i(x_i, u_i), \\ & & \forall i \in \lbrace 0, \ldots, N-1 \rbrace, c_i(x_i, u_i) = 0, \\ & & \forall i \in \lbrace 0, \ldots, N-1 \rbrace, g_i(x_i, u_i) \leq 0, \\ & & c_N(x_N) = 0, \\ & & g_N(x_N) \leq 0 . \end{align}\]

The variables \(x_i, u_i\) represent the states and controls; the functions \(d_i, f_i, c_i, g_i\) are the dynamics, costs, equality constraints, and inequality constraints (respectively); \(s_0\) is the fixed initial state; \(N\) is the number of stages.

1.2 The Regularized Interior Point Method↩︎

The regularized interior point method (IPM) is a way of solving optimization problems of the form

\[\min\limits_{x} f(x) \qquad s.t. \quad c(x) = 0 \wedge g(x) + s = 0 \wedge s \geq 0,\]

where the functions \(f, c, g\) are required to be continuously differentiable.

We start by defining the Barrier-Lagrangian

\[\mathcal{L}(x, s, y, z; \mu) = f(x) - \mu \sum \limits_{i} \log(s_i) + y^T c(x) + z^T (g(x) + s)\]

and the Augmented Barrier-Lagrangian

\[\mathcal{A}(x, s, y, z; \mu, \eta) = \mathcal{L}(x, s, y, z; \mu) + \frac{\eta}{2} (\lVert c(x) \rVert^2 + \lVert g(x) + s\rVert^2).\]

At each iteration of the regularized interior point method, the search direction \((\Delta x, \Delta s, \Delta y, \Delta z)\) is computed by solving the linear system

\[\begin{align} \label{ipm-4x4-newton-kkt} \begin{bmatrix} P & 0 & C^T & G^T \\ 0 & S^{-1} Z & 0 & I \\ C & 0 & -\frac{1}{\eta} I & 0 \\ G & I & 0 & -\frac{1}{\eta} I \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta s \\ \Delta y \\ \Delta z \end{bmatrix} = - \nabla \mathcal{L}(x, s, y, z; \mu), \end{align}\tag{1}\]

where \(P\) is a positive-definite approximation of \(\nabla^2_{xx} \mathcal{L}(x, s, y, z; \mu)\), \(C = J(c)(x)\), \(G = J(g)(x)\), and \(S, Z\) denote the diagonal matrices with entries \(s, z\) respectively.

We call \(x, s\) the primal variables and \(y, z\) the dual variables. As usual, we require that \(s, z > 0\) always holds. Below, we will show that primal search direction determined by the regularized interior point method is guaranteed to be a descent direction of \(\mathcal{A}(\cdot, \cdot, y, z; \mu, \eta)\) at \((x, s)\) unless the primal variables have already converged.

Below, we use \(D( \cdot ; \cdot )\) to represent the directional derivative operator.

Thus, the Augmented Barrier-Lagrangian \(\mathcal{A}(\cdot, \cdot, y, z; \mu, \eta)\) can be used as the merit function for a line search over the primal variables \((x, s)\).

Finally, note that the variable \(\Delta s\) can be eliminated from 1 via \[\Delta s = -Z^{-1} S (\Delta z - \mu S^{-1} e + z) = -Z^{-1} S \Delta z + \mu Z^{-1} e - s,\] where \(e\) represents the all-\(1\) vector, resulting in the linear system

\[\begin{align} \label{ipm-3x3-newton-kkt} \begin{bmatrix} P & C^T & G^T \\ C & -\frac{1}{\eta} I & 0 \\ G & 0 & -(W + \frac{1}{\eta} I) \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \\ \Delta z \end{bmatrix} = - \begin{bmatrix} \nabla_x \mathcal{L}(x, s, y, z; \mu, \eta) \\ c(x) \\ g(x) + \mu Z^{-1} e \end{bmatrix}, \end{align}\tag{2}\] where \(W = Z^{-1} S\).

1.3 Applying the regularized IPM to optimal control↩︎

As shown in 1.2, the line search directions used in the regularized interior point method are computed by solving a linear system of the form

\[\begin{align} \begin{bmatrix} P & C^T & G^T \\ C & -\frac{1}{\eta} I & 0 \\ G & 0 & -W - \frac{1}{\eta} I \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \\ \Delta z \end{bmatrix} = - \begin{bmatrix} r_x \\ r_y \\ r_z \end{bmatrix}, \end{align}\]

where \(\eta > 0\) and \(W\) is a diagonal matrix with positive elements.

Note that \(\Delta z\) can be eliminated via \(\Delta z = (W + \frac{1}{\eta} I)^{-1}(G \Delta x + r_z)\), resulting in

\[\begin{align} \begin{bmatrix} P + G^T (W + \frac{1}{\eta} I)^{-1} G & C^T \\ C & -\frac{1}{\eta} I \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix} = - \begin{bmatrix} r_x + G^T (W + \frac{1}{\eta} I)^{-1} r_z \\ r_y \end{bmatrix}. \end{align}\]

Any component of \(\Delta y\) corresponding to an equality constraint other than the dynamics can be eliminated in the same fashion. Importantly, note that the only constraints that have cross-stage dependencies are the dynamics, so these eliminations preserve the stagewise nature of the problem.

This leaves us with a linear system that is a Dual-Regularized LQR problem.

1.4 Dual-Regularized LQR↩︎

A dual-regularized LQR problem is a linear system of the form \[\begin{align} \label{primal-dual-linear-system} \begin{bmatrix} P & C^T \\ C & -\delta I \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = - \begin{bmatrix} s \\ c \end{bmatrix}, \end{align}\tag{3}\] where \[\begin{align} \label{optimal-control-specialization} P &= \begin{bmatrix} P_0 & & \\ & \ddots & \\ & & P_N \end{bmatrix}, \quad P_i = \begin{cases} \begin{bmatrix} Q_i & M_i \\ M_i^T & R_i \end{bmatrix}, & \text{if } 0 \leq i < N, \\ Q_i, & \text{if } i = N, \end{cases} \\ C &= \begin{bmatrix} -I & & & & & & & & & & \\ A_0 & B_0 & -I & & & & & & & & \\ & & A_1 & B_1 & -I & & & & & & \\ & & & & A_2 & B_2 & -I & & & & \\ & & & & & & & \ddots & -I & & \\ & & & & & & & & A_{N-1} & B_{N-1} & -I \end{bmatrix}, \\ s &= \begin{bmatrix} q_0 \\ r_0 \\ \vdots \\ q_{N-1} \\ r_{N-1} \\ q_N \end{bmatrix}, \quad c = \begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{N-1} \end{bmatrix}, \quad x = \begin{bmatrix} x_0 \\ u_0 \\ \vdots \\ x_{N - 1} \\ u_{N - 1} \\ x_N \end{bmatrix}, \quad y = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_N \end{bmatrix}, \quad \delta > 0. \end{align}\tag{4}\]

Moreover, the matrices \(P_i\) are expected to be symmetric and positive semi-definite, and the matrices \(R_i\) are expected to be symmetric and positive definite.

Note that, when \(\delta = 0\), this linear system is the KKT system associated with a standard LQR problem (i.e. not dual-regularized, and viewed as an optimization problem).

The linear system 3 can be seen as the KKT system associated with the optimization problem \[\begin{align} \label{dual-regularized-lagrangian} \max \limits_{y} \min \limits_{x} \frac{1}{2} x^T P x + s^T x - \frac{\delta}{2} y^T y + y^T (C x + c) . \end{align}\tag{5}\] We call this the dual-regularized Lagrangian. By abuse of notation, when specializing it to 4 , we also call this a dual-regularized LQR problem.

Specializing 5 to the case of discrete-time optimal control problems, as in 4 , this becomes \[\begin{align} \label{min-max-with-ys} \max \limits_{y_0, \ldots, y_N} \min \limits_{x_0, u_0, \ldots, x_N} L(x_0, u_0, \ldots, x_N, y_0, \ldots, y_N), \end{align}\tag{6}\] where \[\begin{align} L(x_0, u_0, \ldots, x_N, y_0, \ldots, y_N) = & \sum \limits_{i=0}^{N-1} \left( \frac{1}{2} x_i^T Q_i x_i + x_i^T M_i u_i + \frac{1}{2} u_i^T R_i u_i + q_i^T x_i + r_i^T u_i \right) + \frac{1}{2} x_N^T Q_N x_N + q_N^T x_N \\ & + y_0^T (c_0 - x_0) + \sum \limits_{i=0}^{N-1} \left( y_{i+1}^T (A_i x_i + B_i u_i + c_{i+1} - x_{i + 1}) \right) - \frac{\delta}{2} \sum \limits_{i=0}^{N} y_i^T y_i . \end{align}\]

1.5 Associative Scans Overview↩︎

Associative scans are a common parallelization mechanism used in functional programming, first introduced in [3]. They were used in [2] to derive a simple method for solving (primal) LQR problems in \(O(\log^2(m) + \log(N) \log^2(n))\) parallel time, where \(N, n, m\) are respectively the number of stages, states, and controls. Note, however, that [2] states that the complexity of their algorithm is \(O(m + \log(N) n)\), apparently stemming from some confusion around the parallel time complexity of matrix inversion (which is \(\log(n)^2\) for \(n \times n\) matrices, as shown in [4]).

Given a set \(\mathcal{X}\), a function \(f: \mathcal{X} \times \mathcal{X} \rightarrow \mathcal{X}\) is said to be associative if \(\forall a, b, c \in \mathcal{X}, f(f(a, b), c) = f(a, f(b, c))\). The forward associative scan operation \(S_f(x_1, \ldots, x_n; f)\) can be inductively defined by \[\begin{align} S_f(x_1; f) &= x_1, \\ S_f(x_1, \ldots, x_{i + 1}; f) &= (y_1, \ldots, y_i, f(y_i, x_{i + 1})), \end{align}\] where \(y_1, \ldots, y_i = S_f(x_1, \ldots, x_i; f)\). Similarly, the reverse associative scan operation \(S_r(x_1, \ldots, x_n; f)\) can be inductively defined by \[\begin{align} S_r(x_1; f) &= x_1, \\ S_r(x_1, \ldots, x_{i + 1}; f) &= (f(x_1, y_2), y_2, \ldots, y_{i + 1}), \end{align}\] where \(y_2, \ldots, y_{i + 1} = S_r(x_2, \ldots, x_{i + 1}; f)\). [3] provides a method for performing associative scans of \(N\) elements in parallel time \(O(\log(N))\).

We will employ associative scans do derive the parallel dual-regularized Riccati recursion, effectively extending [2].

1.6 Motivation↩︎

The reader may ask why we wish to regularize the dynamics and initial state constraints, whose Jacobians always have full row rank due to the presence of the identity blocks. In fact, if we regularized all equality and inequality constraints except for the dynamics, and performed the variable eliminations described in 1.3, we would be left with a standard LQR problem. The issue with doing so is that the computed primal search direction would no longer be guaranteed to be a descent direction of the Augmented Barrier-Lagrangian with the pre-specified penalty parameter.

1.7 Related Work↩︎

The Differentiable Dynamic Programming (DDP) algorithm, introduced in [5], is still the most common method used in practice for numerical optimal control, in spite of the disadvantages associated with single-shooting methods (e.g. the impossibility of independently warm-starting state trajectories).

The Riccati recursion used for solving standard LQR problems was first derived in [1], and was first combined with an interior point method in [6].

More recently, FATROP [7] explores using a different extension of Riccati’s recursion (targeting stagewise equality constraint handling), combining it with an interior-point approach for inequality constraint handling.

An efficient parallel LQR algorithm with logarithmic complexity on the number of stages was first derived in [2], where it was also used to solve unconstrained non-convex discrete-time optimal control problems.

2 Sequential Algorithm↩︎

We start with some simple (but helpful) lemmas that we will employ to eliminate some of the variables in 6 .

For \(0 \in \lbrace 0, \ldots, N \rbrace\), we define the cost-to-go functions \[\begin{align} \label{cost-to-go-value-functions} V_{k}(x_k) = \max \limits_{y_{k+1}, \ldots, y_N} \min \limits_{u_k, x_{k+1}, \ldots, u_{N-1}, x_N} & \sum \limits_{i=k}^{N-1} \left( \frac{1}{2} x_i^T Q_i x_i + x_i^T M_i u_i + \frac{1}{2} u_i^T R_i u_i + q_i^T x_i + r_i^T u_i \right) \\ & + \frac{1}{2} x_N^T Q_N x_n + q_N^T x_n \\ & + \sum \limits_{i=k}^{N-1} \left( y_{i+1}^T \left( A_i x_i + B_i u_i + c_{i + 1} - x_{i + 1} \right) \right) - \frac{\delta}{2} \sum \limits_{i=k}^{N-1} y_{i+1}^T y_{i+1} . \\ \end{align}\tag{7}\]

This brings us to the key theorem.

3 Parallel Algorithm↩︎

First, we can eliminate the control variables \(u_i\) from 6 . Computing the gradients with respect to \(u_i\), equating them to \(0\), and solving for \(u_i\) yields \[\begin{align} u_i = -R_i^{-1} \left( M_i^T x_i + r_i + B_i^T y_{i+1} \right). \end{align}\]

Noting that \[\begin{align} \frac{1}{2} u_i^T R_i u_i &= \frac{1}{2} x_i^T M_i R_i^{-1} M_i^T x_i + r_i^T R_i^{-1} M_i^T x_i + \frac{1}{2} y_{i+1}^T B_i R_i^{-1} B_i^T y_{i+1} + y_{i+1}^T B_i R_i^{-1} M_i^T x_i + r_i^T R_i^{-1} B_i^T y_{i+1} \\ & + \frac{1}{2} r_i^T R_i^{-1} r_i, \\ x_i^T M_i u_i &= -x_i^T M_i R_i^{-1} M_i^T x_i - r_i^T R_i^{-1} M_i^T x_i - y_{i+1}^T B_i R_i^{-1} M_i^T x_i, \\ r_i^T u_i &= -r_i^T R_i^{-1} M_i^T x_i - r_i^T R_i^{-1} r_i - y_{i+1}^T B_i R_i^{-1} r_i, \\ y_{i+1}^T B_i u_i &= -y_{i+1}^T B_i R_i^{-1} B_i^T y_{i+1} - y_{i+1}^T B_i R_i^{-1} M_i^T x_i - r_i^T R_i^{-1} B_i^T y_{i+1}, \end{align}\] and discarding constant terms, the optimization problem 6 becomes \[\begin{align} \label{min-max-without-us} \max \limits_{y_0, \ldots, y_N} \min \limits_{x_0, \ldots, x_N} & \sum \limits_{i=0}^{N-1} \left( \frac{1}{2} x_i^T (Q_i - M_i R_i^{-1} M_i^T) x_i + (q_i - M_i R_i^{-1} r_i)^T x_i \right) \\ & + \frac{1}{2} x_N^T Q_N x_N + q_N^T x_N + \sum \limits_{i=0}^{N-1} \left( y_{i+1}^T \left( \left( A_i - B_i R_i^{-1} M_i^T \right) x_i + \left( c_{i+1} - B_i R_i^{-1} r_i \right) - x_{i + 1} \right) \right) \\ & - \frac{1}{2} \sum \limits_{i=0}^{N-1} y_{i+1}^T \left( \delta I + B_i R_i^{-1} B_i^T \right) y_{i+1} + y_0^T (c_0 - x_0) - \frac{\delta}{2} y_0^T y_0 . \end{align}\tag{20}\]

For \(i \in \lbrace 0, \ldots, N - 1 \rbrace\), we define \[\begin{align} \label{lqr-gpu-init-running} P_{i \rightarrow i+1} &= Q_i - M_i R_i^{-1} M_i^T, \\ p_{i \rightarrow i+1} &= q_i - M_i R_i^{-1} r_i, \\ A_{i \rightarrow i+1} &= A_i - B_i R_i^{-1} M_i^T, \\ C_{i \rightarrow i+1} &= \delta I + B_i R_i^{-1} B_i^T, \\ c_{i \rightarrow i+1} &= c_{i+1} - B_i R_i^{-1} r_i, \\ \end{align}\tag{21}\]

Similarly, we define \[\begin{align} \label{lqr-gpu-init-terminal} P_{N \rightarrow N+1} &= Q_N, \\ p_{N \rightarrow N+1} &= q_N, \\ A_{N \rightarrow N+1} &= 0, \\ C_{N \rightarrow N+1} &= 0, \\ c_{N \rightarrow N+1} &= 0, \\ \end{align}\tag{22}\]

For \(0 \leq i < j \leq N+1\), we define the interval value functions \[\begin{align} V_{i \rightarrow j}(x_i, x_j) = \max \limits_{y_{i+1}, \ldots, y_j} \min \limits_{x_{i+1}, \ldots, x_{j-1}} & \sum \limits_{k=i}^{j-1} \left( \frac{1}{2} x_k^T P_{k \rightarrow k + 1} x_k + p_{k \rightarrow k+1}^T x_k \right) \\ & + \sum \limits_{k=i}^{j-1} \left( y_{k+1}^T \left( A_{k \rightarrow k+1} x_k + c_{k \rightarrow k+1} - x_{k + 1} \right) \right) \\ & - \frac{1}{2} \sum \limits_{k=i}^{j-1} y_{k+1}^T C_{k \rightarrow k+1} y_{k+1} \\ \end{align}\]

Note that, by definition, \[\begin{align} V_i(x_i) = \min \limits_{x_{N+1}} V_{i \rightarrow N+1} (x_i, x_{N+1}). \end{align}\]

Since \[\begin{align} \max \limits_{y_{N+1}} \min \limits_{x_{N+1}} -y_{N+1}^T x_{N+1} = 0, \end{align}\] achieved by forcing \(x_{N+1} = 0\), it follows that \[\begin{align} V_i(x_i) = V_{i \rightarrow N+1} (x_i, 0). \end{align}\]

Thus, 20 amounts to solving \[\begin{align} \max \limits_{y_0} \min \limits_{x_0} V_{0}(x_0) -y_0^T (c_0 - x_0) - \frac{\delta}{2} y_0^T y_0 . \end{align}\]

It was shown in [2] that the functions \(V_{i \rightarrow j}\) admit representations of the form \[\begin{align} V_{i \rightarrow j} (x_i, x_j) = \max \limits_{y} \left( \frac{1}{2} x_i^T P_{i \rightarrow j} x_i + p_{i \rightarrow j}^T x_i - \frac{1}{2} y^T C_{i \rightarrow j} y + y^T \left( A_{i \rightarrow j} x_i + c_{i \rightarrow j} - x_j \right) \right), \end{align}\] modulo constant additive terms that are independent of \(x_i, x_j\) and can thus be discarded.

The following combination rules, established in [2], can be applied to compute \(V_{i \rightarrow k}\) from \(V_{i \rightarrow j}\) and \(V_{j \rightarrow k}\): \[\begin{align} \label{lqr-gpu-combo} P_{i \rightarrow k} &= A_{i \rightarrow j}^T \left( I + P_{j \rightarrow k} C_{i \rightarrow j} \right)^{-1} P_{j \rightarrow k} A_{i \rightarrow j} + P_{i \rightarrow j}, \\ p_{i \rightarrow k} &= A_{i \rightarrow j}^T \left( I + P_{j \rightarrow k} C_{i \rightarrow j} \right)^{-1} \left( p_{j \rightarrow k} + P_{j \rightarrow k} c_{i \rightarrow j} \right) + p_{i \rightarrow j}, \\ A_{i \rightarrow k} &= A_{j \rightarrow k} \left( I + C_{i \rightarrow j} P_{j \rightarrow k} \right)^{-1} A_{i \rightarrow j}, \\ C_{i \rightarrow k} &= A_{j \rightarrow k} \left( I + C_{i \rightarrow j} P_{j \rightarrow k} \right)^{-1} C_{i \rightarrow j} A_{j \rightarrow k}^T + C_{j \rightarrow k}, \\ c_{i \rightarrow k} &= A_{j \rightarrow k} \left( I + C_{i \rightarrow j} P_{j \rightarrow k} \right)^{-1} \left( c_{i \rightarrow j} - C_{i \rightarrow j} p_{j \rightarrow k} \right) + c_{j \rightarrow k}. \\ \end{align}\tag{23}\]

Note that we depart from [2] in 21 , but not in 23 or 22 .

A reverse associative scan [3] can be used to compute the \(P_{i \rightarrow N+1}, p_{i \rightarrow N+1}\) (i.e. the \(P_i, p_i\) in 16 ) in \(O(\log(N) \log (n)^2)\).

Next, the \((I + \delta P_i)^{-1}\) can be computed in \(O(\log(n)^2)\) parallel time (i.e. \(O(1)\) with respect to \(N\)).

We can then compute \(x_0\) in \(O(\log(n))\) parallel time via 19 and the \(K_i, k_i\) from 14 in \(O(\log(n) + \log (m)^2)\) parallel time (i.e. \(O(1)\) with respect to \(N\) in both cases).

Finally, we wish to compute the \(u_i, x_{i+1}\) from the \(K_i, k_i\) via 11 and 14 . Note that the sequential LQR forward pass has \(O(N)\) parallel time complexity. However, as done in [2], we can reduce the computation of the \(x_i\) to a sequential composition of affine functions, which can also be parallelized with an associative scan [3] in \(O(\log(N) \log(n))\) parallel time. This will be described in 3.1. Due to 2,

\[\begin{align} x_{i+1} & = (I + \delta P_i)^{-1} \left( A_i x_i + B_i u_i + c_{i+1} - \delta p_{i+1} \right) \\ &= (I + \delta P_i)^{-1} \left( A_i x_i + B_i \left( K_i x_i + k_i \right) + c_{i+1} - \delta p_{i+1} \right) \\ &= (I + \delta P_i)^{-1} \left( A_i + B_i K_i \right) x_i + (I + \delta P_i)^{-1} ( B_i k_i + c_{i+1} - \delta p_{i+1}) . \end{align}\]

Once these affine functions have been composed, they can be independently applied to \(x_0\) to recover all the \(x_i\) in \(O(\log(n))\) parallel time. The \(u_i\) can then be computed in \(O(\log (\max (m, n)))\) parallel time by independently evaluating \(u_i = K_i x_i + k_i\). Note that both of these parallel times are \(O(1)\) with respect to \(N\).

Thus, the combined parallel time complexity of our method is \(O(\log(m)^2 + \log(N) \log(n)^2)\).

3.1 Composing Affine Functions with Associative Scans↩︎

In this section, we describe an algorithm for composing \(N\) affine functions \(F_i(x) = M_i x_i + m_i\), as done in [2], with \(O(\log(N) \log(n))\) parallel time.

Letting \(\mathcal{X} = \mathbb{R}^{n} \times \mathbb{R}^{n \times n}\) and letting \(f: \mathcal{X} \rightarrow \mathcal{X}\) be defined by \(f((a, B), (c, D)) = (Da + c, DB)\), we claim that \(f\) is associative. This is simple to verify: \[\begin{align} &f(f((a, B), (c, D)), (e, F)) \\ = &f((Da + c, DB), (e, F)) \\ = &(F(Da + c) + e, F(DB)) \\ = &((FD)a + Fc + e, (FD)B) \\ = &f((a, B), (Fc + e, FD)) \\ = &f((a, B), f((c, D), (e, F))). \end{align}\]

Moreover, note that \(f\) is the affine function composition operator, as \(D(Bx + a) + c = (DB)x + (Da + c)\). Thus, it suffices to apply a forward associative scan to \(f\) to obtain all compositions \(F_0, F_1 \circ F_0, \ldots, F_{N-1} \circ \cdots \circ F_0\) in \(O(\log(N) \log(n))\) parallel time.

4 Software Contributions↩︎

We provide JAX implementations of both the sequential and parallel dual-regularized LQR algorithms [8], including unit tests verifying the correctness of the solutions provided by our methods on a set of random examples satisfying the required definiteness properties.

Modern interior-point solvers, such as IPOPT [9], while supporting several generic linear system solvers, are not designed to easily accommodate linear solvers that specialize to specific problem types (e.g. optimal control problems). In fact, until recently, IPOPT installations did not include all of the required headers for users to integrate any custom user-side linear system solvers, showing that this was not being done at all.

In order to address this limitation, we released a simple regularized interior point solver in [10], which allows users to easily plug in callbacks specializing certain operations to their specific problem type, including:

• a KKT system factorization callback;

• a KKT system solve callback;

• a KKT system residual computation callback.

This effectively splits the solver into a shared backend and a problem-specific frontend. We do exactly that for the case of optimal control in [11].

Adding to this, we also implemented an integration with QDLDL [12] supporting arbitrary user-provided KKT permutations (for optimal fill-in prevention) in [13], resulting in an efficient sparse interior point solver that avoids dense operations entirely. Alternatively, a sparse linear algebra code generation library, such as SLACG [14], can be used to solve the KKT systems and compute their residuals; this often achieves substantial speed-ups for solving sparse problems, as no index computations are performed at runtime.

Examples can be found in [15]. In all cases, with correct usage, dynamic memory allocation is entirely avoided.

All of the provided packages are free and open-source (MIT-licensed).

5 Conclusion↩︎

We derived extensions of the sequential and parallel Riccati Recursions for dual-regularized LQR problems, allowing us to easily and efficiently numerically solve constrained non-linear discrete-time optimal control problems via the regularized interior point method. We also provided free, efficient, open-source implementations of the described algorithms, which are fully unit-tested on random well-defined examples, establishing the correctness of our derivations and implementations.

References↩︎