Efficient Online Update of Model Predictive Control in Embedded Systems Using First-Order Methods

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ Esta es la versión aceptada del artículo publicado en: This is a accepted manuscript of a paper published in: IEEE Control Systems Letters (2023): 11 diciembre 2023 DOI: 10.1109/LCSYS.2023.3341283 Copyright: ©2023IEEE El acceso a la versión publicada del artículo puede requerir la suscripción de la revista. Access to the published version may require subscription. “© 20XX 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” This article has been accepted for publication in IEEE Control Systems Letters. This is the author’s version which has not been fully edited and content may change prior to final publication. Citation information: DOI https://doi.org/10.1109/LCSYS.2023.3341283 1 Efficient online update of model predictive control in embedded systems using first-order methods Victor Gracia, Pablo Krupa, Teodoro Alamo, Daniel Limon Abstract—Model Predictive Control (MPC) is typically characterized for being computationally demanding, as it requires solving optimization problems online; a particularly relevant point when considering its implementation in embedded systems. To reduce the computational burden of the optimization algorithm, most solvers perform as many offline operations as possible, typically performing the computation and factorization of its expensive matrices offline and then storing them in the embedded system. This improves the efficiency of the solver, with the disadvantage that online changes on some of the ingredients of the MPC formulation require performing these expensive computations online. This article presents an efficient algorithm for the factorization of the key matrix used in several firstorder optimization methods applied to linear MPC formulations, allowing its prediction model and cost function matrices to be updated online at the expense of a small computational cost. We show results comparing the proposed approach with other solvers from the literature applied to a linear time-varying system. I. INTRODUCTION Model Predictive Control (MPC) [1], [2] is a well-known control policy due to its ability to optimize the system operation while satisfying its constraints. However, its control law is computationally demanding, as it requires solving an optimization problem every sample time. Many linear MPC solvers have been proposed in the literature using different optimization methods, such as active set [3], interior point [4], or various first-order methods (FOM) such as gradient descent [5], Douglas-Rachford [6], the alternating direction method of multipliers (ADMM) [7] or FISTA [8], to name a few. In particular, FOMs are a popular choice for their use in embedded systems [9], [10], due to their simplicity, good practical performance, and the availability of certification results on their iteration complexity [11], [12], [13]. Due to the low-resource nature of embedded systems, the implementation of these methods often avoids any non-strictly necessary online computations by computing all possible solver components offline, which are then stored in the embedded system [12], [14], [15]. This approach is particularly advantageous when the computation of these elements requires computationally expensive operations, such as matrix inversion or factorization. In particular, factorization of the key solver ingredients is an approach typically taken to reduce the computational complexity per-iteration of the algorithm, see, e.g., [4], [12], [14], [16]. The disadvantage, however, is that in a practical setting it may Department of Systems Engineering and Automation, Universidad de Sevilla, Spain (e-mails: [email protected],[email protected],[email protected], [email protected]). Corresponding author: Victor Gracia. The authors acknowledge support from Grants PID2022-141159OB-I00 and PID2022-142946NA-I00, funded by MCIN/AEI/ 10.13039/501100011033 and by “ERDF A way of making Europe”, and Grant PDC2021-121120-C21 funded by MCIN/AEI/10.13039/501100011033 and by the “European Union NextGenerationEU/PRTR”. no longer be possible to change online some of the ingredients of the MPC controller, such as its prediction model or cost function matrices. That is, the operations required to update the solver for the new resulting optimization problem may be too computationally expensive for the embedded system given the sample time, mostly due to the previously mentioned matrix computations and factorization. A typical scenario that requires online adjustment of some of these elements is when dealing with linear time-varying (LTV) or linear parameter-varying (LPV) systems [17], where the prediction model of the MPC controller is typically updated every sample time. In fact, the need to update the prediction model will be a requirement when implementing any adaptive MPC control scheme [18]. Additionally, being able to tune online the cost function matrices of the MPC controller is a welcome feature when working in a practical setting, either to fine tune the controller online to improve its performance when controlling the real system or to adapt to a changing economic criterion on the “best” way to control it. Some quadratic programming (QP) solvers provide the option to recompute and refactor its computationally expensive ingredients online, e.g., [14], [15]. However, the required computation time can be significant when compared to their solve time. In [19], a Ricatti recursion is used to compute the factorization of the linear system to be solved at each iteration of the algorithm. The computation time of this factorization procedure does not depend on the prediction horizon of the MPC controller if its terminal cost is taken as the solution of the discrete-time Ricatti equation associated to the infinitehorizon unconstrained problem. In this article, we present an algorithm for updating the ingredients used in several FOM-based solvers, including the popular methods ADMM and FISTA, when applied to standard linear MPC formulations whose prediction model and/or weight matrices change online.This algorithm, in essence, is an efficient computation of the Cholesky decomposition of the main computationally expensive matrix used in several solvers [12], [20]. The key feature of the proposed approach is that the matrix factorization is performed exploiting the particular structure that arises from many standard linear MPC formulations, providing a computational advantage when compared to the online factorization routines provided by solvers for generic QPs from the literature. This may lead to tractable refactorization of the solver ingredients in cases where the refactorization schemes of generic QP solvers are too slow due to the short sample time or the limitations of the embedded system where the solver is implemented. To illustrate this, we show numerical results comparing the proposed approach to several state-of-the-art QP solvers applied to an LTV system. ©2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. arXiv:2309.07996v2 [eess.SY] 16 Dec 2023 This article has been accepted for publication in IEEE Control Systems Letters. This is the author’s version which has not been fully edited and content may change prior to final publication. Citation information: DOI https://doi.org/10.1109/LCSYS.2023.3341283 2 The article is structured as follows. Section II presents the problem formulation. In Section III, we present the proposed algorithm for the online computation and factorization of the solver ingredient. Section IV shows the numerical results. Finally, Section Vsummarizes the main conclusions. Notation: The symbol Zb adenotes the set of integers from the integers ato b, both included. For a matrix A, its element in row iand column jis denoted as Ai,j . We denote ∥x∥Q. =px⊤Qx. The symbol Indenotes the identity matrix of dimension n. Given two vectors xand y,x≤(≥)y denotes componentwise inequalities. For vectors x1to xN, (x1, x2, . . . , xN)denotes the column vector formed by their concatenation. Given scalars and/or matrices A1, A2, . . . , AN (not necessarily of the same dimensions), we denote by diag(A1, A2, . . . , AN)the block diagonal matrix formed by the concatenation of A1to AN. II. PROBLEM FORMULATION Consider a discrete-time system described by x(k+ 1) = Ax(k) + Bu(k),(1) where x(k)∈Rnand u(k)∈Rmare the state and input at sample time k, respectively, subject to box constraints x(k)≤x(k)≤x(k),(2a) u(k)≤u(k)≤u(k),(2b) where the bounds x(k),x(k)∈Rnand u(k),u(k)∈Rm satisfy x(k)< x(k)and u(k)< u(k). The control objective is to steer the system to a reference pair (xr, ur), assumed to be an admissible equilibrium point of the system. This objective can be attained using one of many linear MPC formulations proposed in the literature [1], [2], [12]. For instance, consider the standard linear MPC formulation min x,u N−1 X j=0 ||xj−xr||2 Q+||uj−ur||2 R+||xN−xr||2 T(3a) s.t. xj+1 =Axj+Buj, j ∈ZN−1 0,(3b) x0=x(k),(3c) xj≤xj≤xj, j ∈ZN 1,(3d) uj≤uj≤uj, j ∈ZN−1 0,(3e) where Nis the prediction horizon; Q∈Rn×n,R∈Rm×m and T∈Rn×nare positive definite matrices; and the decision variables xand ucontain the state and input variables, respectively, along the prediction horizon. It is well known that, under mild assumptions, system (1) under the control law u(k) = u∗ 0, where u∗ 0is the optimal solution of the decision variable u0of (3), asymptotically converges to the admissible reference (xr, ur)without violating the system constraints [2]. Therefore, the implementation of the MPC controller requires obtaining the optimal solution of (3) at each sample time. We remark that formulation (3) has been selected as an example instead of any other linear MPC formulation in order to provide the reader with concrete expressions. However, as shown later on, the proposed update procedure can be adapted to other linear MPC formulations with minor modifications. Problem (3) can be rewritten as a standard QP problem min z 1 2z⊤Hz +q⊤z(4a) s.t. Gz =b, (4b) z≤z≤z, (4c) by taking H=diag(R, Q, R, . . . , Q, R, T),(5a) q=−(Rur, Qxr, Rur, . . . , Qxr, Rur, Txr),(5b) G=       B−In0· · · · · · 0 0A B −In· · · 0 0 0 .........0 0 0 0 A B −In       ,(5c) b=−(Ax(k),0,0, ..., 0),(5d) z= (u0, x1, u1, ..., xN−1, uN−1, xN),(5e) z= (u0, x1, u1, ..., xN−1, uN−1, xN).(5f) There are many optimization methods that can efficiently solve problem (4) (active-set method, interior point methods, FOMs, etc.). Among them, FOMs have received a significant amount of attention from the control community, especially when dealing with the implementation of MPC in embedded systems, whose restricted computational and memory resources are a good match for these simple-to-implement algorithms with low resource requirements. Furthermore, a major advantage of solving (4) using FOMs is that the ingredients of the optimization method typically have simple sparse structures due to the sparse, banded nature of (5a) and (5c); a fact exploited by several QP solvers from the literature to obtain efficient implementations, e.g., [12], [14], [16]. For instance, in [12] the authors propose sparse solvers based on ADMM or FISTA for problem (4) arising from simple standard MPC formulations such as (3). The proposed solvers take advantage of the particular structures of Hand Gthat arise from these MPC formulations, i.e., (5a) and (5c). Specifically, they require solving a linear system Wz =dat each iteration of the algorithm, where W. =G(H+ρI)−1G⊤ for some non-negative scalar ρwhose value depends on the specific FOM being used. The solvers in [12] solve this linear system by backward-forward substitution using the Cholesky decomposition of W, i.e., the matrix Wcsatisfying W=W⊤ cWc, whose simple banded structure is exploited to provide an efficient algorithm. In particular, the structure of Wcresulting from problem (4) is given by Wc=          β1α10· · · 0 0β2α2· · · 0 . . .. . .β3.... . . . . .. . .. . ....αN−1 0 0 0 · · · βN          ,(6) where βk∈Rn×n,k∈ZN 1, are upper triangular matrices and αk∈Rn×n,k∈ZN−1 1, are dense. Matrices αkand βk, which ©2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. This article has been accepted for publication in IEEE Control Systems Letters. This is the author’s version which has not been fully edited and content may change prior to final publication. Citation information: DOI https://doi.org/10.1109/LCSYS.2023.3341283 3 Algorithm 1: Computation of the sub-matrices of Wc Require: Igredients that affect W(A,B,Q,R,T) 1Compute the matrices in Γ(8) 2for k= 1 to N−1do 3Compute βkusing (10a) and (10b) 4Compute αkusing (10c) 5end for 6Compute βNusing (10a) and (10b) Output: β1to βN, α1to αN−1 depend on the values of the MPC ingredients A,B,Q,Rand T, are the main computationally-expensive ingredients of the solver, and are thus typically computed offline and stored in the embedded system, where they are used to efficiently solve the system of equations Wz =d[21, §5.1]. A limitation of these FOM solvers, therefore, is that changing some of the MPC ingredients online requires recomputing and/or refactoring certain matrices, which can be rather timeconsuming. In the previous example, changing any of the ingredients A,B,Q,Ror Twould require an online computation of Wc(6). In other solvers from the literature, such as [14], [15], [16], similar sparse matrices need to be recomputed and refactored if similar MPC ingredients are changed. The following section presents an efficient algorithm for updating the βkand αksub-matrices of Wc(6) obtained from QP problems (4) whose Hand Gingredients have banded structures similar to the ones shown in (5a) and (5c). The motivation behind this algorithm is to provide optimization solvers which can efficiently deal with MPC ingredients that change online, such as when dealing with LTV or LPV systems instead of the time-invariant system (1). We note that MPC formulation (3) is taken as an example whose associated Wc matrix has the structure shown in (6), but that most other standard linear MPC formulations lead to nearly identical structures, requiring, as we will show, only minor modifications of the procedure presented in the following section. III. ONLINE COMPUTATION OF THE SOLVER INGREDIENTS We now present equations for recursively computing submatrices αkand βkof Wc(6) associated with the QP problem (4) for the MPC formulation (3). That is, we provide expressions for computing the non-zero sub-matrices of the Cholesky decomposition of matrix W. =G(H+ρI)−1G⊤. We present the final expressions of the equations. The appendix provides an intuitive overview of how they were obtained. Let us introduce the notation ˆ Q. = (Q+ρI)−1,ˆ R. = (R+ρI)−1,ˆ T. = (T+ρI)−1,(7a) Y. =Aˆ QA⊤, Z . =Bˆ RB⊤,(7b) XQ. =Y+Z+ˆ Q, XT. =Y+Z+ˆ T, (7c) where we recall that the value of ρ≥0depends on the FOM being used [12]. Setting α0= 0 for convenience, we denote Γk. =     Z+ˆ Q, if k= 1, XQ,if k∈ZN−1 2, XT,if k=N, (8a) (8b) (8c) γk i,j . = i−1 X l=1 βk l,iβk l,j + n X q=1 αk−1 q,i αk−1 q,j .(9) In (8), Γkcan be seen intuitively as a variable that gathers the elements of the diagonal band of Wthat affect βk. In (9), γk groups terms of the non-diagonal elements of βkand elements of αk−1that affect subsequent elements of βk. Keeping this in mind, the nonzero elements βk i,j of βk∈Rn×n, k ∈ZN 1, and the elements αk i,j of αk∈Rn×n, k ∈ZN−1 1, are given by βk i,i =qΓk i,i −γk i,i, i ∈Zn 1,(10a) βk i,j =Γk i,j −γk i,j βk i,i , j ∈Zn i+1,(10b) (βk)⊤αk=−(Aˆ Q)⊤.(10c) Note that the expressions in (7a) require matrix inversions, which are cheap when Q, R and Tare small and/or diagonal. However, if this is not the case, one could directly provide ˆ Q, ˆ Rand ˆ Tin order to reduce its computation time. The computation of αkusing (10c) is simple, since each of its columns is computed as the solution of the upper-triangular linear system (βk)⊤αk [j]=−(Aˆ Q)⊤ [j], where we use subindex [j]to indicate column jof the matrix. Algorithm 1shows the procedure for computing αkand βk. The computational complexity of Algorithm 1when applied to the Wresulting from (3) is O(m3+n3N)flops. Remark 1. The computational complexity of Algorithm 1can be reduced if certain additional assumptions are made, such as matrices Qand Rbeing diagonal. An interesting point about equations (9) and (10) is that they are a consequence of the structure of matrix Wcshown in (6). That is, equations (9) and (10) are valid for any MPC formulation whose solver requires solving a linear system Wz =dwhose Cholesky factorization Wchas the structure shown in (6). Thus, the only difference between Algorithm 1 applied to each formulation is the matrices required to determine the value of Γk(8) and, possibly, the matrix on the right-hand side of (10c). For the MPC formulation (3), the values of Γkfor each kare shown in (8) and (7), but for other MPC formulations these expressions may differ. Indeed, let us now consider the standard MPC formulation with terminal equality constraint [12, Eq. (8)], given by min x,u N−1 X j=0 ||xj−xr||2 Q+||uj−ur||2 R(11a) s.t. xj+1 =Axj+Buj, j ∈ZN−1 1,(11b) x0=x(k),(11c) xj≤xj≤xj, j ∈ZN 1,(11d) uj≤uj≤uj, j ∈ZN−1 0,(11e) xN=xr.(11f) In this case, matrices Hand Gof the resulting QP problem (4) are nearly identical to the ones shown in (5a) and (5c); the only differences being in the final rows of said matrices ©2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. This article has been accepted for publication in IEEE Control Systems Letters. This is the author’s version which has not been fully edited and content may change prior to final publication. Citation information: DOI https://doi.org/10.1109/LCSYS.2023.3341283 4 due to its different terminal ingredients. Thus, equations (8)- (10) are still valid for updating its αkand βkmatrices, with the exception that (8c) now uses Y+Zinstead of XT. This change is a consequence of the inclusion of the terminal equality constraint and the removal of the terminal cost in (11), which results in a different expression for ΓN. However, the expressions for the other submatrices of Wcdo not change. For other standard linear MPC formulations whose Hand Gmatrices maintain banded structures similar to the ones shown in (5a) and (5c), similar small modifications of Γwill be required, but overall Algorithm 1will remain the same. This is the case, for instance, with the MPC formulations from [4, §V.A], [10], [12], [19]. Furthermore, we find that this same procedure may also apply to some non-standard MPC formulations, such as in the MPC for tracking solver presented in [20] (c.f., [20, Eq. (11)]). Remark 2. The expression of ΓNfor (11)differs from the one for (3)because (11f)adds an equality constraint, affecting the structure of (5a)and (5c), and thus W. However, if we consider the more general terminal constraint xN∈ Xf, where Xfis a polyhedron, inequality constraints are added instead of an inequality constraint. Therefore, Algorithm 1would still be valid as long as the inequality constraints are handled so as to not affect the banded structure of Wexploited by Algorithm 1, e.g., uncoupling the equality and inequality constraints using ADMM as in [12] or by using a suitable dualization approach of the inequality constraints as in [19]. Remark 3. An alternative approach to the one taken here would be to implement a generic sparse Cholesky decomposition procedure for sparse matrices. The use of generic sparse matrix factorization routines is used by some QP solvers to provide an efficient way to update the solver ingredients online, such as in the OSQP solver [14], which implements a generic sparse LDL factorization routine. This approach, while efficient due to its sparse implementation, is for generic sparse matrices, i.e., using generic sparse matrix representations such as the “compressed sparse column”. In contrast, our approach does not make use of sparse matrix representations, but instead directly computes the non-zero elements of the Cholesky decomposition of Wby exploiting the specific structure of the matrices Hand Gshown in (5a)and (5c). That is, our approach loses generality in favor of improving its efficiency by being particularized to the specific structure that typically arises in standard linear MPC formulations, thus providing better computation times than the ones obtained using the update procedures of generic QP solvers. Remark 4. The expressions provided in (7)-(10)for the computation of αk, k ∈ZN−1 1, and βk, k ∈ZN 1, can be easily adapted in order to handle variations in A,B,Qand R between prediction steps, i.e., having Aj,Bj,Qj,Rj,j∈ZN 1, in problem (3). In this case, (7)includes the inversion of N matrices, i.e., computing ˆ Q−1 j,ˆ R−1 j,j∈ZN 1, as well as the matrix multiplications required to compute Yj=Ajˆ QjA⊤ j and Zj=Bjˆ RjB⊤ j, j ∈ZN 1. Consequently, the number of operations of Algorithm 1increases. Particularly, it presents a computational complexity of O((m3+n3)N)flops. IV. NUMERICAL RESULTS We consider the continuous stirred-tank reactor (CSTR) system from [22], where a reactant Ais transformed into a product Bby an exothermic and irreversible chemical reaction. The objective is to control the reaction volume (V), the concentration of A(CA), and the temperature inside the reactor (T). The manipulable variables are the output flow rate (qs) and the input coolant flow rate (qc). The non-linear continuous-time equations of the system are given by dV (t) dt =qe−qs(t), dCA(t) dt =qe V(t)(CAe −CA(t)) −k0e −E RT (t)CA(t), dT(t) dt =qe V(t)(Te−T(t)) −k1e −E T(t)CA(t) +qc(t) V(t)k2(1 −e −k3 qc(t))(Tce −T(t)), where k1=∆Hk0 ρCp,k2=ρCpc ρcCp,k3=hA ρcCpc and the values of the system parameters are presented in [22]. We consider an LPV model of the system, obtained from the above nonlinear equations as shown in [22], with a frozen scheduling parameter, i.e., the matrices of the discrete LPV system model change at every sample time, and thus also the MPC prediction model, which remains constant along its prediction horizon. We take the sampling time of the system as 0.1seconds, and control it using the MPC formulation (3). We implement the ADMM and FISTA solvers from [12], available in the SPCIES toolbox (version v0.3.8) [23], which include Algorithm 1. We also solve the MPC optimization problem using the OSQP (version 0.6.2) [14] and SCS (version 3.2.3) [15] solvers, which include a generic update procedure, as discussed in Remark 3. We note that OSQP also implements ADMM and SCS implements a similar operator splitting algorithm. The simulations are run on an Intel Core i5-1135G7. Table Ishows the computation times when taking average results from 200 simulations using random initial conditions and references within the system constraints described in [22]. Computation times are separated into “update-time”, which is the time taken by the procedure to update the solver ingredients due to the modifications of the MPC parameters and ingredients (this includes the computation time of Algorithm 1 in the case of the SPCIES solvers), and “solve-time”, which is the remaining computation time. The results show a clear difference in terms of computation time between the updatetime of Algorithm 1w.r.t. OSQP or SCS. In order to show that the updating procedure we propose works properly when combined with time-varying MPC approaches, and that the time results we shown in Table Iare obtained from a suitably designed MPC controller (in the sense of proper closed-loop behavior), we include Figure 1, which shows the closed-loop result of one of the simulations of Table Iusing the solver from [23], where the MPC model changes every sample time due to the LPV approach used, and the cost function matrices Qand Rwere changed from Q= 10Into Q=Inand from R=diag(10,0.1) to R=Im ©2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. This article has been accepted for publication in IEEE Control Systems Letters. This is the author’s version which has not been fully edited and content may change prior to final publication. Citation information: DOI https://doi.org/10.1109/LCSYS.2023.3341283 5 Update-time [microseconds] Solve-time [microseconds] Solver NAverage Median Maximum Minimum Average Median Maximum Minimum % FISTA 5 1.03 (0.08) 1.02 (0.07) 1.30 (0.17) 0.97 (0.06) 2.69 2.06 12.25 1.36 27.75 ADMM 5 1.05 (0.11) 1.04 (0.10) 1.36 (0.20) 0.96 (0.08) 71.01 23.49 686.3 15.14 1.45 OSQP 5 18.03 18.01 19.81 17.45 20.41 18.39 43.29 17.59 46.91 SCS 5 101.6 101.6 109.9 99.20 180.1 212.2 245.7 104.1 36.08 FISTA 15 1.68 (0.20) 1.67 (0.19) 2.11 (0.57) 1.62 (0.17) 6.14 4.47 29.40 2.71 21.52 ADMM 15 1.81 (0.22) 1.80 (0.21) 2.35 (0.47) 1.72 (0.19) 156.5 100.2 1428.3 50.83 1.15 OSQP 15 36.95 36.93 38.73 35.94 53.39 47.54 159.65 46.05 40.90 SCS 15 233.7 233.7 244.9 229.6 560.8 597.7 663.6 339.1 29.42 The update times of algorithms FISTA and ADMM have been obtained using SPCIES with Algorithm 1in its update procedure. In brackets, we include the update time of the same SPCIES solvers when Algorithm 1is not executed, i.e., without updating A,B,Qand/or R. Column “%” shows (Avg. update-time)/(Avg. update-time + Avg. solve-time). The ADMM penalty parameter of SPCIES is set to ρ= 0.01. TABLE I: Computation times of the MPC solvers applied to the CSTR system for the MPC formulation (3). T Set-point (a) State T. qc Set-point Bounds (b) Input qc. Fig. 1: Non-linear CSTR model controlled with LPV frozen scheduling parameter MPC approach. SPCIES OSQP SCS (a) Dependency on n(m= 2, N = 5) (b) Dependency on m(n= 60, N = 5) (c) Dependency on N(n= 4, m = 2) Fig. 2: Dependency of update times on n,mand N. at time t= 9 seconds, explaining the different system behavior between the first and second half of the simulation Finally, Figure 2presents results on the computation time of Algorithm 1in terms of the number of states (n), inputs (m) and the prediction horizon (N) as well as the update times when using the OSQP and SCS solvers. Note that Figure 2c is shown in logarithmic scale in order to improve its readability. The results show that the proposed approach provides faster computation times than when using state-of-the-art generic QP solvers, as well as better scaling as the dimensions of the MPC problem increase. This is due to our particularization of the update procedure to the specific matrix structure that arises in linear MPC problems. We note that the update-time of Algorithm 1grows linearly with the prediction horizon N, as discussed in Section III, since the number of matrices αk and βkis linear in N. Remark 5. We do not include a comparison with [19] because we focus on comparisons with state-of-the-art solvers available online. The computational complexity of the “updatestep” [19, Algorithm 3] does not depend on the prediction horizon Nif the terminal cost of the MPC controller is taken as the solution of the discrete LQR Riccati equation. In this case, its computation time would be smaller than the one obtained with Algorithm 1. However, if this is not the case, or if the system is not LTI along the prediction horizon, Algorithm 1 outperforms it in terms of the number of required operations. V. CONCLUSION This article has presented an efficient method to perform the matrix factorization required by several FOM optimization solvers from the literature for linear MPC formulations. The proposed update method is based on exploiting the particular structure that typically arises from standard MPC formulations, thus providing a computational advantage w.r.t. the general sparse update procedures used by many state-of-the-art QP solvers, as shown in the numerical results we have presented, where computation times were up to an order of magnitude faster. Combined with the solvers presented in [12] (available ©2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. This article has been accepted for publication in IEEE Control Systems Letters. This is the author’s version which has not been fully edited and content may change prior to final publication. Citation information: DOI https://doi.org/10.1109/LCSYS.2023.3341283 6 in [23]), which are also based on exploiting MPC-specific matrix structures, this becomes an interesting option for applying MPC in under-powered devices when an online change of the prediction model or the tuning parameters of the MPC controller are required. APPENDIX This section provides an intuitive overview of how the update expressions for αkand βkpresented in Section III are derived. Recalling the definition (7), W. =G(H+Iρ)−1G⊤ is given by W=         Z+ˆ Q−ˆ QA⊤0· · · 0 −Aˆ Q XQ−ˆ QA⊤0 0 0.........0 0 0 −Aˆ Q XQ−ˆ QA⊤ 0 0 0 −Aˆ Q XT         .(13) The structure of its Cholesky decomposition is given by (6). Using the identity W=W⊤ cWc, we can derive the following equations relating αkand βkwith the matrices in (13): (β1)⊤β1=Z+ˆ Q, (β1)⊤α1=−(Aˆ Q)⊤, (β2)⊤β2=XQ−(α1)⊤α1,(β2)⊤α2=−(Aˆ Q)⊤, . . . (βN)⊤βN=XT−(αN−1)⊤αN−1. For matrices βk, we can use their upper-triangular nature to further expand the above equations to obtain expressions for their non-zero elements. For instance, for β1, we have (β1)⊤β1=     β1 1,10 0 . . ....0 β1 1,n · · · β1 n,n         β1 1,1· · · β1 1,n 0.... . . 0 0 β1 n,n     =Z+ˆ Q, from where we derive (β1 1,1)2= (Z+ˆ Q)1,1, β1 1,2= (Z+ˆ Q)1,2/β1 1,1, . . . , β1 1,n = (Z+ˆ Q)1,n/β1 1,1,(β1 2,2)2= (Z+ˆ Q)2,2−(β1 1,2)2, β1 2,3= [(Z+ˆ Q)2,3−β1 1,2β1 1,3]/β1 2,2, . . . , β1 2,n = [(Z+ˆ Q)2,n −β1 1,2β1 1,n]/β1 2,2, . . . , (β1 n,n)2= (Z+ˆ Q)n,n −(β1 1,n)2− · · · − (β1 n−1,n)2. (14) Continuing with β2, we have (β2)⊤β2=XQ−(α1)⊤α1, which leads to the same expressions as (14) but substituting β1 with β2and Z+ˆ Qwith XQ−(α1)⊤α1. The same procedure can be performed for each βkup to βN, which satisfies (βN)⊤βN=XT−(αN−1)⊤αN−1. Equations (10a) and (10b) are a generalization of (14), where Γk(8) accounts for the matrices of (13) that affect each βk, e.g., Z+ˆ Qfor β1or XTfor βN, and γ(9) accounts for the effect of the “previous” values of the elements of βkand αk−1that affect the current βk i,j, as shown in (14). Finally, we take the convention α0= 0 so that we can also add α0to the expression for β1in (10). REFERENCES [1] E. F. Camacho and C. B. Alba, Model Predictive Control. Springer Science & Business Media, 2013. [2] J. B. Rawlings, D. Q. Mayne, and M. Diehl, Model predictive control: theory, computation, and design, 2nd ed. Madison, Wisconsin: Nob Hill Publishing, 2017. [3] D. Arnstr¨ om, A. Bemporad, and D. Axehill, “A dual active-set solver for embedded quadratic programming using recursive LDLTupdates,” IEEE Transactions on Automatic Control, vol. 67, no. 8, pp. 4362–4369, 2022. [4] A. Domahidi, A. U. Zgraggen, M. N. Zeilinger, M. Morari, and C. N. Jones, “Efficient interior point methods for multistage problems arising in receding horizon control,” in 51st IEEE Conference on Decision and Control (CDC), 2012, pp. 668–674. [5] T. Skibik and M. M. Nicotra, “Analysis of time-distributed model predictive control when using a regularized primal–dual gradient optimizer,” IEEE Control Systems Letters, vol. 7, pp. 235–240, 2023. [6] R. Halvgaard, L. Vandenberghe, N. K. Poulsen, H. Madsen, and J. B. Jørgensen, “Distributed model predictive control for smart energy systems,” IEEE Transactions on Smart Grid, vol. 7, no. 3, pp. 1675–1682, 2016. [7] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends® in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011. [8] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 1, pp. 183–202, 2009. [9] P. Giselsson, “Improved fast dual gradient methods for embedded model predictive control,” IFAC Proceedings Volumes, vol. 47, no. 3, pp. 2303– 2309, 2014. [10] D. Kouzoupis, H. J. Ferreau, H. Peyrl, and M. Diehl, “First-order methods in embedded nonlinear model predictive control,” in European Control Conference (ECC), 2015, pp. 2617–2622. [11] P. Zometa, M. K¨ ogel, T. Faulwasser, and R. Findeisen, “Implementation aspects of model predictive control for embedded systems,” in American Control Conference (ACC), 2012, pp. 1205–1210. [12] P. Krupa, D. Limon, and T. Alamo, “Implementation of model predictive control in programmable logic controllers,” IEEE Transactions on Control Systems Technology, vol. 29, no. 3, pp. 1117–1130, 2021. [13] P. Krupa, O. Inverso, M. Tribastone, and A. Bemporad, “Certification of the proximal gradient method under fixed-point arithmetic for boxconstrained QP problems,” Automatica, vol. 159, 2024. [14] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “OSQP: An operator splitting solver for quadratic programs,” Mathematical Programming Computation, vol. 12, no. 4, pp. 637–672, 2020. [15] B. O’Donoghue, “Operator splitting for a homogeneous embedding of the linear complementarity problem,” SIAM Journal on Optimization, vol. 31, pp. 1999–2023, August 2021. [16] M. Garstka, M. Cannon, and P. Goulart, “COSMO: A conic operator splitting method for convex conic problems,” Journal of Optimization Theory and Applications, vol. 190, no. 3, pp. 779–810, 2021. [17] J. Bokor and G. Balas, “Linear parameter varying systems: A geometric theory and applications,” IFAC Proceedings Volumes, vol. 38, no. 1, pp. 12–22, 2005, 16th IFAC World Congress. [18] G. Tao, “Multivariable adaptive control: A survey,” Automatica, vol. 50, no. 11, pp. 2737–2764, 2014. [19] P. Patrinos and A. Bemporad, “An accelerated dual gradient-projection algorithm for embedded linear model predictive control,” IEEE Transactions on Automatic Control, vol. 59, no. 1, pp. 18–33, 2014. [20] P. Krupa, I. Alvarado, D. Limon, and T. Alamo, “Implementation of model predictive control for tracking in embedded systems using a sparse extended ADMM algorithm,” IEEE Transactions on Control Systems Technology, vol. 30, no. 4, pp. 1798–1805, 2022. [21] P. Krupa, “Implementation of MPC in embedded systems using first order methods,” Ph.D. dissertation, University of Seville, 2021, available at arXiv:2109.02140. [22] E. Bernardi, H. A. Pipino, and E. J. Adam, “Full-order output observer applied to a linear parameter varying system with unknown input,” in IEEE Congreso Bienal de Argentina (ARGENCON), 2020, pp. 1–8. [23] P. Krupa, V. Gracia, D. Limon, and T. Alamo, “Spcies: Suite of predictive controllers for industrial embedded systems,” 2020. [Online]. Available: https://github.com/GepocUS/Spcies ©2023 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information.