{'fraction_non_alphanumeric': 0.04119629978049545, 'fraction_numerical': 0.01724678582627783, 'mean_word_length': 4.245271381578948, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Low-Rank Modifications of Riccati Factorizations for Model Predictive Control

In Model Predictive Control (MPC) the control input is computed by solving a constrained finite-time optimal control (CFTOC) problem at each sample in the control loop. The main computational effort is often spent on computing the search directions, which in MPC corresponds to solving unconstrained finite-time optimal control (UFTOC) problems. This is commonly performed using Riccati recursions or generic sparsity exploiting algorithms. In this work the focus is efficient search direction computations for active-set (AS) type methods. The system of equations to be solved at each AS iteration is changed only by a low-rank modification of the previous one, and exploiting this structured change is important for the performance of AS type solvers. In this paper, theory for how to exploit these low-rank changes by modifying the Riccati factorization between AS iterations in a structured way is presented. A numerical evaluation of the proposed algorithm shows that the computation time can be significantly reduced by modifying, instead of re-computing, the Riccati factorization. This speed-up can be important for AS type solvers used for linear, nonlinear and hybrid MPC.


Introduction
Model Predictive Control (mpc) is a control strategy where the applied control input is computed by minimizing a cost function while satisfying constraints on the states and control inputs. It has become one of the most widely used advanced control strategies in industry, and some important reasons for its popularity are that it can handle multivariable systems and constraints on states and control inputs in a structured way, [1]. In each sample of the mpc control loop a constrained finite-time optimal control (cftoc) problem is solved on-line, which creates a need for efficient optimization algorithms. Note that similar linear algebra is also useful in off-line applications such as explicit mpc solvers. The mpc problem and the corresponding cftoc problem can be of different types depending on which system and problem formulation that is used. Some common variants are linear mpc, nonlinear mpc and hybrid mpc. In many cases the main computational effort when solving the cftoc problem boils down to compute the search directions, which corresponds to solving unconstrained finite-time optimal control (uftoc) problems. The uftoc problems can be solved using for example Riccati recursions, and some examples of how optimization routines have been sped up by using Riccati recursions are [2,3,4,5,6,7,8,9,10,11,12,13].
The use of Riccati recursions in active-set (as) methods for optimal control can be found as early as in [2]. In this reference a Riccati recursion is used to factor the major block of the kkt matrix, and for the other block standard lowrank modifications of factorizations are used on a dense system of equations of the size of the number of active inequality constraints. The computational complexity of this algorithm grows quadratically in the number of active inequality constraints. An alternative sparse non-Riccati factorization is used in [14], and the factorization is updated after changes in the as iterations.
In as methods it is often crucial to modify the factorization of the kkt matrix instead of re-factorizing it between as iterations, [15]. Since this has traditionally not been considered possible when using the Riccati factorization, it has sometimes been argued that this factorization is not suitable for as methods. However, in [11] a method for making low-rank modifications of the Riccati factorization by exploiting the structured changes between as iterations was introduced, showing that this is indeed possible. The work in [11] is limited to problems with non-singular control input weight matrices and simple control input bounds, and modifications of the kkt matrix is only possible at a single time index.
The main contribution in this paper is the extension of the theory in [11] to handle more general forms of uftoc problems, where the kkt matrix can be singular. The derivation of this result looks similar to the one in [11], but here more technical depth is added since additional mathematical tools are needed in this paper to cope with the singularity of the kkt matrix. In this paper it is also described how to modify the factorization for more general modifications of the kkt matrix where constraints are simultaneously added (or removed) at different time indices. Both these generalizations can be important when using for example dual projection as solvers like the one in [8]. Furthermore, in [11] only bound constraints on the control inputs are considered in the cftoc problem, whereas it will be shown in this paper how the theory can be applied to problems with both state and control input constraints. A more detailed description of the material presented in this paper can be found in the thesis in [16].
In this article, S n ++ (S n + ) denotes symmetric positive (semi) definite matrices with n columns, Z i,j = {z ∈ Z | i ≤ z ≤ j}, and R (A) denotes the range space of a matrix A.
Section 2 introduces the cftoc problem and Section 3 optimization preliminaries. In Section 4 the main result is presented, and in Section 5 it is shown how to use this for more general problems. The numerical results and conclusion are presented in Section 6 and Section 7, respectively.

Problem Formulation
For linear mpc problems, the corresponding cftoc problem consists of a quadratic objective function and affine dynamics constraints. For now, consider only upper and lower bounds on the control signal. Let t denote the time-index in the mpc optimization problem (i.e., t = 0 is the current time), N the prediction horizon, x t ∈ R nx the state, u t ∈ R nu,t the control input,x ∈ R nx the initial state, and the stacked states and control inputs, respectively. The cftoc problem can then be written in the form min. x,u where Q x,N ∈ S nx + and Q t ∈ S nx+nu,t + . Let the matrix Q t and the vector l t ∈ R nx+nu,t be partitioned as Note that the more common additional assumption Q u,t ∈ S nu,t ++ is not used in this problem formulation in order to, for example, be able to represent dual mpc problems.
Furthermore, define λ t+1 as the dual variable for the equality constraint −x t+1 + A t x t + B t u t + a t = 0 and µ t as the dual variable for the inequality constraint

Optimization Preliminaries
The cftoc problem (2) is a convex quadratic programming (qp) problem. Hence, it can be solved using several different types of optimization methods, where one common type is as methods, see for example, [17,15].

Active-set QP methods
as methods solve a qp problem by determining the set of constraints that are active, i.e., hold with equality, at the optimal solution. This set of active constraints is denoted the optimal active set, and an as solver operates by finding this set of constraints iteratively, [15]. Since the optimal active set is usually not known a priori, an as solver starts with an initial set of constraints that are forced to hold with equality, and then changes this so-called working set by adding or removing constraints until the optimal active set has been determined. These modifications of the working set are usually relatively small and the modifications of the corresponding kkt matrix between as iterations are thus of low rank. The modification techniques presented in this work can be used both by traditional as solvers where one constraint is added or removed to the working set at each iteration, and for solvers that add or remove several constraints to the working set at each iteration such as those presented in [15,8,9]. Let W j denote the subset of the working set that contains the indices of the inequality constraints that temporarily hold with equality at as iteration j, and let W c j denote the set of inequality constraints that are temporarily disregarded at as iteration j. In problem (2) only control input constraints are used, and hence forcing a constraint to hold with equality corresponds to removing that control input as an optimization variable from the optimization problem and treating it as a constant. Similarly, by disregarding an inequality constraint, the corresponding control input becomes unconstrained. This can be formalized by introducing w t as the free part of the control inputs and v t as the fixed part as follows where Π t is a permutation matrix satisfying Π T t Π t = I. Here u t (W j ) is used to denote the control inputs at time t with corresponding constraints in W j . Using this notation, B t , Q u,t , Q xu,t andl u,t can be partitioned consistently with w t and v t : By using this partitioning of the control input and the corresponding matrices, the uftoc problem that is solved at as iteration j to compute the search direction is given by min. x,w where Computing the sequence of search directions in an as type solver applied to the cftoc problem (2) hence corresponds to solving a sequence of uftoc problems in the form (7).

Standard Riccati recursion
The solution to the uftoc problem (7) is computed by solving a set of linear equations known as the kkt optimality conditions. The special structure of the uftoc problem considered in this work corresponds to a sparse, almost block diagonal, kkt system which can be solved very efficiently using a Riccati recursion, see, e.g., [2,3,8,16]. The Riccati recursion consists of a factorization of the kkt matrix (Algorithm 1), followed by back-and forward substitutions (algorithms 2-4) for computing the solution to (7), [5]. Algorithm 1, which is the computationally demanding part of the Riccati recursion, computes the variables P t+1 ∈ S nx + and K t+1 ∈ R nw,t×nx using the auxiliary variables where M t+1 ∈ S and H t+1 ∈ R nx×nw,t by construction. Since Q w,t ∈ S nw,t + it follows that also G t+1 ∈ S nw,t + . When one (or more) G t+1 is singular a non-unique Riccati factorization still exists, but the solution of the kkt system is either non-unique or non-existing, [8,16]. How to handle this case is determined at the solver level, and one way to do this is presented in [8].
Using the Riccati recursion to compute the solution to (7) requires O (N ) complexity, compared to O N 3 or O N 2 for dense solvers that re-factorize, or modify the factorizations of, the kkt matrix without exploiting the uftoc problem structure, respectively. For more information, see, e.g., [15].

Low-Rank Modification of the Riccati Factorization
A standard approach to improve the performance of an as solver is to modify the factorization of the kkt matrix instead of re-factorizing it between as iterations, [15]. Here it will be shown how to modify the Riccati factorization (Algorithm 1) between as iterations when solving a cftoc problem (2). Since Q u,t ∈ S nu,t + , and hence possibly also Q w,t ∈ S nw,t + , the kkt matrix for the uftoc problem that is solved to compute the search direction can be singular (some G t+1 in Algorithm 1 can be singular). The derivations in this section are similar to the one in [11], but the extension presented here adds more technical

6:
Compute and store a factorization of G t+1 .

7:
Compute a solution K t+1 to: Compute a solution k t+1 to: depth since additional mathematical tools such as generalized Schur complements (gscs), the quotient formula for gscs and the Moore-Penrose pseudoinverse are required to cope with the possibly singular kkt matrix. For a detailed description of these, see for example [18,19,20]. Furthermore, in [11] it was only shown how to modify the Riccati factorization when modifying the working set at a single time instance. If constraints at different time indices are added or removed in the same as iteration, the factorization can be modified by performing a sequence of complete modifications. However, in this section it will be shown how to handle either adding or removing several constraints at different time indices by instead gradually increasing the size of the modification of the factorization. If the solver both adds and removes constraints in the same as iteration, the factorization must be modified sequentially. Note that as the size of the modification increases, it might be better to re-compute the remaining part of the factorization from scratch. Which approach that is most efficient depends on for example the size of the modification, and can be investigated off-line. That work is however outside the scope of this paper.
By introducing the gsc operator as / † , the gsc with respect to G t+1 of Here G † t+1 is the Moore-Penrose pseudo inverse of G t+1 . Hence, by using Line 7 in Algorithm 1 and basic properties of the pseudo inverse, the matrix P t in Line 8 in Algorithm 1 can be calculated as Lemma 1 (Quotient formula for gsc). Let the positive semi-definite matrices M 0 andM 0 be partitioned as Then Proof. Lemma 1 follows directly from Theorem 4 in [19]. The details are given in the proof of Lemma 4.1 in [16].
In this paper, a tilde will be used to indicate a matrix that has been modified. Hence, the modified version of some matrix X is denotedX.

Sequence of low-rank modifications
Assume thatP t+1 for some t ∈ Z 0,N −1 is a downdate of P t+1 , given by (the superscript "−" indicates a downdate) with Later in this section, and in sections 4.2 and 4.3, Lemma 1 will be used to show that this assumption holds for all modifications presented in this paper. The downdate is considered to be of low rank ifk < n x . It will now be shown how this downdate of P t+1 affects the matrices in the Riccati factorization for the time-steps τ ∈ Z 0,t . By substituting P t+1 in lines 3-5 in Algorithm 1 withP t+1 from (14), straightforward calculations giveF The equations in (15) can be written in matrix form as SinceP t+1 ∈ S nx + andM t+ is defined as in (9), it follows thatM t+ is positive semidefinite by construction. Now, definê and letM t+1 be the second diagonal block ofM t+1 . Note thatG t+1 =M t+1 / † C − t+1 . From [18] it follows thatM t+1 0 since C − t+1 0 and V − t+1 T ∈ R C − t+1 by assumption, andM t+1 =M t+1 / † C − t+1 0 in (16). Hence alsoM t+1 0, and Lemma 1 thus can be used to computeP t =M t+1 / †G t+1 =M t+1 / †M t+1 (first equality in (12)). By using the block partitioning , the modified version of P t is computed using the second equality in (12) in Lemma 1 as where Using similar calculations, an update of P t+1 in the form with C + t+1 ∈ Sk + , V + t+1 ∈ R nx×k and V + t+1 T ∈ R C + t+1 , can be shown to result in the updateP with (hereG t+1 andH t+1 are defined similarly as in (15)) Note that the modified matricesH t+1 andG t+1 are used in (22). Hence, a modification of P t+1 of at most rankk results in a similar modification of P t of (also) at most rankk.
Theorem 1. Consider a modification of at most rankk of P tm ∈ S nx + in Algorithm 1 at a single time instant t m ∈ Z 1,N in either of the forms where C − tm , C + tm ∈ Sk + , V − tm , V + tm ∈ R nx×k , and V − tm T ∈ R C − tm , V + tm T ∈ R C + tm , respectively. Then it holds for all t ∈ Z 0,tm−1 that P t ∈ S nx + is modified as with C − Proof. From the derivations of (18) and (21) it follows that a modification of at most rankk of P t+1 at an arbitrary t ∈ Z 0,N −1 in either of the forms (14) or (20) results in a similar modification of at most rankk of P t . Since P tm is modified as in (23), the proof follows by induction.
The modifiedK t+1 can be computed by solving For the common case whereG t+1 ∈ S nw,t ++ , it is possible to use the Sherman-Morrison-Woodbury formula for efficient computations. For the details the reader is referred to, e.g., [21]. Remark 1. The factorization ofG t+1 is computed by modifying the factorization of G t+1 . Hence, a solution to (25) can be computed without having to re-factorizeG t+1 , which requires less computations than re-solving (25) from scratch.

Removing control input constraints from the working set
Assume that P t+1 is modified as in (14) with a modification of dimensionk. Furthermore, assume that k control input constraints that are affecting the control input at time t are removed from the working set, i.e., temporarily disregarding these constraints that previously were forced to hold. This affects the uftoc problem (7) in the same way as adding k new control inputs. Note that this combination of modifications is more general than the one used in [11], where only modifications of the working set at a single time index was considered. Assume without loss of generality that the new control inputs are appended at the end of w t . Then the matrices B w,t , Q w,t and Q xw,t are modified as givingB w,t ∈ R nx×(nw,t+k) ,Q w,t ∈ S nw,t+k + andQ xw,t ∈ R nx×(nw,t+k) . From lines 4-5 in Algorithm 1 it follows thatG t+1 ∈ S nw,t+k + andH t+1 ∈ R nx×(nw,t+k) are given byG where In analogy with Section 4.1,P t can be computed asP t =M t+1 / †G t+1 , wherẽ M t+1 is computed as in (16) but instead usingG t+1 andH t+1 from (27). By definingM withM t+1 as the second diagonal block, it is clear thatM Hence, as in Section 4.1, Lemma 1 states thatP t = M t+1 / †M t+1 . Using the partitioning A = F t+1 and D = G t+1 (B, C, E and F consistently) ofM t+1 , the second equality in Lemma 1 gives where (31b) Hence, removing k control input constraints at time t from the working set when a modification in the form (14) of P t+1 is already present results in a modification of P t in the same form as (14) but of increased dimensionk + k.

Remark 2.
Note that ifk + k is close to, or larger than, n x it might be better to re-compute the factorization. This trade-off can be investigated off-line by modifying the factorization for different sizes of modifications and determine which alternative is faster, but the details are left as future work.
Remark 4. For the common case whenG t+1 ∈ S nw,t+k ++ , low-rank modifications can be exploited by using the Sherman-Morrison-Woodbury formula for efficient computations. The factorization ofG t+1 is modified as is mentioned in Remark 1.
When removing k constraints from the working set also k components of v t are removed. Hence, also straightforward modifications of B v,t , Q xv,t , Q v,t , Q wv,t ,l v,t are made. However, these changes do not affect the matrices in the factorization and are not presented here, see [16].

Adding control input constraints to the working set
Assume that P t+1 is modified as in the form (20), and that k control input constraints that are affecting the control input at time t are added to the working set at AS iteration j. Adding constraints corresponds to removing these control inputs from the problem and treating them as constants. The impact from this modification on P t is similar to when constraints are removed. Assume, without loss of generality, that the k control inputs are removed from the k last entries of w t . The modified matricesB w,t ,Q w,t andQ xw,t are then obtained from The implicit relations betweenF t+1 ,G t+1 ,H t+1 , F t+1 , G t+1 and H t+1 are therefore given by g,g 0 andh are computed from g, g 0 and h in G t+1 and H t+1 . Note that the modified matrices are on the right hand side.
Here,M t+1 andM t+1 are defined analogously as in (29), but using the matrices in (33). Hence, from Lemma 1 where Hence, adding k control input constraints at time t to the working set when a modification in the form (20) is already present results in a modification of P t in the same form as (20) but of increased dimensionk + k. The modifiedK t+1 can be computed by solving (25), but using the modified matrices in (33). Note that Remark 2 and 4 apply here as well.

Remark 5.
If there is no modification of P t+1 , then C +

Algorithms for modifying the Riccati factorization
Let t m be the largest time index where W j is modified. The theory presented in this section is summarized in Algorithm 5, which starts by modifying the matrices in the factorization according to Section 4.2 or Section 4.3 depending on whether constraints are removed or added to the working set, respectively. Since P tm+1 is not modified Remark 3 or Remark 5, respectively, applies. For t < t m the matrices in the factorization are modified as in Section 4.1, 4.2 or 4.3 depending on the type of modification at time t. Note that only adding or removing constraints is possible at the same as iteration. As is mentioned in Remark 1 standard methods for modifying the factorization of G t+1 should be used to avoid re-computing the factorization. See for example [21,22,15] for details on techniques for modifying factorizations.
For an example with non-singular Q u,t where k constraints are removed at time t m and Cholesky factorizations of G t+1 are used, the computational complexity when modifying the Riccati factorization instead of re-computing it is reduced from approximately O N (n 3 w + n 3 x + n 2 w n x + n 2 x n w ) to approximately O t m (n 2 w n x + n 2 w + kn w n x + kn 2 x ) . If the solution to (25) is computed using the Sherman-Morrison-Woodbury formula the complexity is further reduced to approximately O t m (n 2 w + kn 2 w + kn w n x + kn 2 x ) . Note that the complexity is now linear in t m and quadratic in n x and n w , which shows the gains of modifying the Riccati factorization instead of re-computing it. However, the exact expression for the complexity depends on for example the choice of factorization and modification techniques in Algorithm 1 and 5. To compute the solution to the modified uftoc problem, the recursions in algorithms 2-4 need to be re-computed. Since the factorization is only modified for t ∈ Z 0,tm the backward recursion in Algorithm 2 only needs to be re-computed for t ≤ t m using the modified matricesB w,t ,G t ,H t andP t .

Extension to General Constraints
The cftoc problem arising in many mpc problems in industry often includes constraints on the states and the possibility to control only certain states or a linear combination of states [1], and is of a more general type of problem than (2). Here it will be described how the theory presented in Section 4 can be used to compute the search directions even when solving more general cftoc problems than (2) using an as type solver. Note that the purpose with this section is not to present a complete as solver, but to explain how the theory can be used when solving more general problems than (2). In this section, the superscripts "p" and "d" denote variables related to the primal problem and the dual problem, respectively. and the quadratic terms in the objective function satisfy All computations were performed on an Intel Xeon W3565 @3.2 GHz processor running Linux (version 2.6.32-504.12.2.el6.x86 64) and Matlab (version 9.1.0.441655 (R2016b)). The default settings have been used for the code generation in Matlab, with the exception that the compilation flag '-O3' has been used to optimize the code for speed.
The algorithms are compared by solving random uftoc problems in the form (7), where n x and n w are logarithmically spaced in Z 10,200 . The computation times are averaged over 20 different problems of the same dimensions. In Fig. 1 the computation times are normalized w.r.t. to the maximum computation time 0.65 seconds for the standard Riccati recursion. Here, a problem with N = 10 has been solved after removing a constraint at either t m = 0 (modifying one step of the factorization) or t m = N − 1 (modifying the full factorization), which are the best and worst case for the modifying algorithm, respectively. Furthermore, in Fig. 2 the performance gains for different N and t m ∈ {1, N 4 , N 2 , 3N 4 , N } are investigated by plotting the ratio between the computation times when modifying the Riccati factorization and re-computing it for problems of dimension n x , n w = 10, n x , n w = 103 and n x , n w = 200, respectively. From the figures it is clear that modifying the Riccati factorization instead of re-computing it can significantly reduce the computation time for solving the uftoc problem (7), especially for large problem sizes and/or when only a small part of the factorization is modified. The O (t m ) complexity (independent on N ) result in Section 4.4 is numerically verified in Fig. 2, where it is shown that the performance is similar for N ∈ {10, 20, 40, 60, 80, 100}. The accuracies of the numerical solutions have been measured as the Euclidean norm of the kkt residual for the uftoc problem (7). For a problem with N = 100 and n x , n w = 200 the maximum residual norm is in the order 10 −10 for both the standard Riccati recursion and the proposed algorithm.

Conclusions
This work presents theory and algorithms for modifying the Riccati factorization instead of re-computing it after low-rank modifications of the kkt system have been made. This is possible by exploiting the special structure from the mpc problem, and it can be used to significantly improve the performance of as type solvers by modifying the Riccati factorization between as iterations instead of re-computing it. The algorithm has been evaluated using a C implementation generated from Matlab's code generation framework, and it is shown that significant gains in terms of performance can be obtained using the proposed algorithm. The result shows that Riccati recursions can be employed in as methods without sacrificing the important possibility to exploit low-rank modifications of the kkt systems when computing the search directions required to solve a cftoc problem.