Implementation of MPC in embedded systems using first order methods
Full text
Escuela Superior T´ecnica de Ingenier´ıa Departamento de Ingenier´ıa de Sistemas y Autom´atica Doctoral Thesis Implementation of MPC in embedded systems using first order methods Pablo Krupa Garc´ıa Supervised by: Daniel Limon Marruedo Teodoro Alamo Cantarero Seville, June 2021
Contents Acknowledgements i Notation, conventions and definitions iii 1 Introduction 1 1.1 Motivation and objectives . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Outline of the dissertation . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Publications............................... 5 I Restart schemes for accelerated first order methods 7 2 Preliminaries: Accelerated first order methods 9 2.1 First order methods . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.1 Proximal gradient method . . . . . . . . . . . . . . . . . . . 10 2.1.2 Alternating direction method of multipliers . . . . . . . . . 16 2.1.3 The extended alternating direction method of multipliers . 17 2.2 Accelerated first order methods . . . . . . . . . . . . . . . . . . . . 18 2.2.1 Fast proximal gradient method (FISTA) . . . . . . . . . . . 18 2.2.2 A monotone variant of FISTA . . . . . . . . . . . . . . . . . 23 2.3 Conclusions............................... 24 3 Restart schemes for accelerated first order methods 25 3.1 A brief review of the literature . . . . . . . . . . . . . . . . . . . . 27 3.1.1 Fixed-rate restart schemes . . . . . . . . . . . . . . . . . . . 28 3.1.2 Adaptive restart schemes . . . . . . . . . . . . . . . . . . . 30 3.2 Implementable restart schemes with linear convergence . . . . . . . 34 3.2.1 A restart scheme for FISTA . . . . . . . . . . . . . . . . . . 34 3.2.2 A gradient based restart scheme for FISTA . . . . . . . . . 42 3.2.3 Restart scheme for accelerated first order methods . . . . . 47 3.3 Numericalresults............................ 59 3.3.1 Application to Lasso problems . . . . . . . . . . . . . . . . 60 3.3.2 Application to QP problems . . . . . . . . . . . . . . . . . . 63 3.4 Conclusions and future lines of work . . . . . . . . . . . . . . . . . 66 II Implementation of MPC in embedded systems 67 4 Preliminaries 69 4.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.2 A brief introduction to model predictive control . . . . . . . . . . . 71 4.3 A brief state of the art . . . . . . . . . . . . . . . . . . . . . . . . . 73
5 Sparse solvers for model predictive control 77 5.1 Various structure-exploiting solvers . . . . . . . . . . . . . . . . . . 78 5.1.1 Solving systems of equations with banded decomposition . 78 5.1.2 Solving equality-constrained QPs with banded structure . . 80 5.1.3 Solving box-constrained separable QPs . . . . . . . . . . . . 81 5.2 Solving QPs with FISTA through duality . . . . . . . . . . . . . . 82 5.2.1 QP problem’s dual formulation . . . . . . . . . . . . . . . . 83 5.2.2 Solving the QP’s dual problem with FISTA . . . . . . . . . 85 5.3 Solving QPs with ADMM . . . . . . . . . . . . . . . . . . . . . . . 87 5.4 Simple standard MPC formulations . . . . . . . . . . . . . . . . . . 88 5.4.1 FISTA-based solver for standard MPC . . . . . . . . . . . . 91 5.4.2 ADMM-based solver for standard MPC . . . . . . . . . . . 92 5.5 MPC with terminal quadratic constraint . . . . . . . . . . . . . . . 94 5.5.1 ADMM solver for MPC with terminal quadratic constraint 95 5.5.2 Computation of admissible ellipsoidal invariant sets . . . . 100 5.6 MPCfortracking............................102 5.6.1 Recasting the MPCT formulation for EADMM . . . . . . . 104 5.6.2 EADMM-based solver for the MPCT formulation . . . . . . 105 5.7 TestBenches ..............................108 5.7.1 Chemical plant: double reactor and separator . . . . . . . . 108 5.7.2 Ballandplate..........................111 5.7.3 Oscillating masses . . . . . . . . . . . . . . . . . . . . . . . 112 5.8 Numericalresults............................113 5.8.1 Comparison between the proposed ADMM-based solvers . . 114 5.8.2 Comparison between the ADMMand FISTA-based solvers 118 5.8.3 Standard MPC subject to terminal equality constraint . . . 121 5.8.4 Standard MPC without terminal constraint . . . . . . . . . 122 5.8.5 MPC subject to terminal quadratic constraint . . . . . . . . 123 5.8.6 MPC for tracking . . . . . . . . . . . . . . . . . . . . . . . . 125 5.8.7 Restart methods applied to the FISTA-based solvers . . . . 128 5.9 Conclusions and future lines of work . . . . . . . . . . . . . . . . . 135 6 Harmonic based model predictive control for tracking 137 6.1 The harmonic MPC formulation . . . . . . . . . . . . . . . . . . . 138 6.2 Recursive feasibility of the HMPC formulation . . . . . . . . . . . 144 6.3 Asymptotic stability of the HMPC formulation . . . . . . . . . . . 147 6.4 Selection of the base frequency . . . . . . . . . . . . . . . . . . . . 153 6.5 Advantages of the HMPC formulation . . . . . . . . . . . . . . . . 155 6.5.1 Performance advantages of the HMPC formulation . . . . . 155 6.5.2 Domain of attraction . . . . . . . . . . . . . . . . . . . . . . 159 6.6 Conclusions and future lines of work . . . . . . . . . . . . . . . . . 162 Bibliography 163 List of Algorithms 177 Index 179
Acknowledgements First and foremost, I want to thank my supervisors, Professors Daniel Lim´on and Teodoro ´ Alamo. The passion and energy of Daniel, along with his goodhearted disposition and ever-present willingness to discuss my problems and questions, hooked me in and kept me going. It was in his office, in 2015, that this journey began when I met him to ask about a graduate thesis he was offering. I found myself, a solid hour and a half later, having received a private lecture on MPC simply because I showed some interest. Looking back at that meeting, I have come to realize that it was one of the most life-changing moments of my life. I wouldn’t be here without the support and the knowledge he has shared with me since then. Teodoro’s attitude towards science and research, his strive for excellence, and his seemingly endless knowledge have inspired me and made me grow throughout these years. His constant search of a better algorithm and a better proof, and his obsession with using the proper punctuation marks after every single equation, have, for better or worse, been deeply drilled into my mind. I will always fondly remember the hours spent in his office discussing optimization algorithms and MPC; oftentimes with me struggling to keep up with his line of thought, but finally managing to do so through his sheer perseverance. I don’t know if I would have the patience. I am deeply in debt with both of them; more, I think, than they are aware of. I would also like to thank all the other members of the GEPOC research group; those currently there and those now working elsewhere. They have always been helpful and kind. One cannot ask for a better work environment. I would also like to acknowledge Nilay Saraf, who is a pleasure to work with. I know few people as willing as he to work the extra hour to get a better scientific result. A huge thanks to my family, friends and loved ones. I will not list them all for fear leaving someone out, but there are a few I feel deserve a special mention. To my parents and brother, for the numerous times they have asked about my progress during these years, even though they may not understand half of what I’m talking about. Likewise, I would like to thank my friends Manuel and Isa, who are always supportive (and as insistent as humanly possible). In particular, I would like to thank Manuel for his huge effort in trying to understand my work and all the interest he has shown over the years. Thanks to Laura, who has also always shown a great interest in what I do; going so far as to promise that she will read this dissertation. I’m holding you to your word. A deep and special thanks to Paula, for always being there and encouraging me, particularly during these last few months while I was writing this dissertation. Finally, thanks to all the people I worked with in the Mitsubishi Electric Research Laboratories. In particular, a special thanks to Dr. Claus Danielson, who welcomed and helped me during my stay and from whom I learnt a lot. Pablo Krupa Garc´ıa, Seville, June 2021. i
ii
Notation, conventions and definitions We list here the notation, basic definitions, conventions and well-known mathematical background that we use throughout the paper. Additionally, sectionspecific definitions and notation may be presented throughout the manuscript when needed. The definitions and conventions shown here are standard in the literature, and therefore no proofs nor specific references are provided. We refer the reader to [1, 2, 3] as general references containing most of the definitions and conventions stated here. Spaces and general set notation Let Cand Dbe two sets. We use the standard set notations C ⊂ D (Cis a strict subset of D), D ⊆ D (Cis a subset or equal to D), C ∩D (intersection of Cand D) and C∪D (union of Cand D). If xis an element of set C, we write x∈ C. The empty set is denoted by ∅. The set of real numbers is denoted by R. The set of extended real numbers, i.e., Rextended with −∞ and +∞, is denoted by R. We denote by R≥0and R>0 the set of non-negative and (strictly) positive real numbers, respectively. We use analogous notations for the non-positive (R≤0) and negative (R<0) real numbers as well as for their extended real number counterparts (R>0,R≥0,R<0,R≤0). We denote by [a, b] the set of real numbers xsatisfying a≤x≤b. We use a rounded bracket (instead of a square one), to indicate strict inequality, e.g., (a, b] denotes the set of real numbers xsatisfying a<x≤b. The set of integer numbers is denoted by Z. For integers i, j ∈Zsatisfying i≤j, we denote Zj i . ={x∈Z:i≤x≤j}. We use the same notation used with real numbers to indicate the sets of positive (Z>0), non-negative (Z≥0), negative (Z<0) and non-positive (Z≤0) integer numbers. Vectors and matrices We denote by Rnthe set of n-dimensional real vectors. For x∈Rn, we denote by x(i)its ith component. The inner product (or dot product) of two vectors x, y ∈Rnis denoted by hx, yi=Pn i=1 x(i)y(i). All vectors are considered column vectors unless otherwise specified. For vectors x1∈Rn1, x2∈Rn2, . . . , xN∈RnN, with N∈Z>0, we denote by (x1, x2, . . . , xN) the column vector formed by the concatenation of column vectors x1to xN. Vector (x1, x2, . . . , xN)∈Rn1×Rn2× ···×RnNcan also be viewed as the Cartesian product of vectors x1to xN. The transposed of a vector x∈Rnis denoted by x>. For vectors x, y ∈Rn, the notations x≤y,x>y, etc. are to be taken componentwise. We denote by Rm×nthe space of mby ndimensional real matrices. For M∈ Rm×n, we denote by M(i,j)its (i, j)th component, M>its transposed and M−1its iii
inverse (if Mis non-singular). We say that Mis diagonal if M(i,j)= 0 whenever i6=j. The space of diagonal matrices in Rn×nis denoted by Dn. We denote by Sn ++,Sn +⊂Rn×nthe spaces of positive definite and positive semi-definite matrices in Rn×n, respectively. We denote their diagonal counterparts by Dn ++ ⊂Sn ++ and Dn +⊂Sn +. Given a symmetric matrix M, its maximum and minimum eigenvalues are given by λmax(M) and λmin(M), respectively. For matrices and/or scalars M1, M2, . . . , MN(not necessarily of the same dimension), with N∈Z>0, we denote by diag(M1, M2, . . . , MN) the diagonal concatenation of M1to MN. We denote by 1n×m∈Rn×mand 0n×m∈Rn×mthe matrices of all ones and all zeros, respectively. The vectors of ones and zeros in Rnare denoted by 1nand 0n, respectively. The identity matrix of dimension nis denoted by In. We may drop the dimensions of the matrices/vectors if they are clear from the context, and instead simply write 1,0or I. Norms Given a∈R,|a|denotes its absolute value. Given x∈Rn, we denote: the standard Euclidean norm by kxk2. =phx, xi; the M-weighted Euclidean norm by kxkM. =phx, Mxi(provided that M∈Sn ++); the `1-norm by kxk1. =Pn i=1 |x(i)|; the `∞-norm by kxk∞. = maxi∈Zn 1|x(i)|; and kxkto represent any vector norm. The spectral norm of a matrix M∈Rn×mis denoted by kMk. =pλmax(M>M). Definition N.1 (Dual norm).Given any norm k·k :Rn→R≥0, its dual norm k·k∗:Rn→R≥0is defined as kxk∗. = supz∈Rn{hx, zi:kzk ≤ 1}. If k·k is an M-weighted Euclidean norm kxk=phx, Mxi, where M∈Sn ++, then its dual norm is given by kxk∗=phx, M−1xi. Sequences A sequence of elements xjindexed by j∈ J ⊆ Zis denoted by {xj}, where the set Jwill be clear from the context. Additionally, we use {xj}≥0if J=Z≥0, and {xj}>0if J=Z>0. If the index set Jis finite, then we denote by a bold xthe Cartesian product of {xj}. We often times append the index term jas a superscript, instead of a subscript, i.e., we write xjto denote element jof an ordered sequence. In particular, we do this to express the iterates of some algorithms. For signals (or sequences) that evolve over time we use the following specific notation. Let x∈Rnbe a vector whose value changes over time, which we index by letter t. We denote by x(t) the value of xat time instant t. iv
Set definitions Definition N.2 (Convex set).A set C ⊆ Rnis convex if αx + (1 −α)y∈ C, ∀x, y ∈ C,∀α∈[0,1]. Definition N.3 (Affine set).A set C ⊆ Rnis said to be affine if it contains all the lines that pass through pairs of points x, y ∈ C with x6=y. Definition N.4 (Closed, open, bounded and compact sets).A non-empty set C ∈ Rnis closed if the limit of every converging sequence {xk}, with xk∈ C, is also contained in C. It is open if its complement {x∈Rn:x6∈ C} is closed. It is bounded is there exists c∈R>0such that kxk ≤ c,∀x∈ C. It is compact if it is closed and bounded. Definition N.5 (Affine hull).The affine hull of a set C, which is denoted by aff(C), is the intersection of all the affine sets containing C. Definition N.6 (Interior and relative interior of a set).A point x∈Rnis an interior point of a set C ⊆ Rnis there exists an open sphere S. ={z∈Rn: kz−xk2< }that is contained in C. The set of all interior points of Cis called the interior, and is denoted by int(C). A point x∈Rnis a relative interior point of a set C ⊆ Rnis there exists an open sphere Sfor some ∈R>0such that S∩aff(C)⊂ C. The set of all relative interior points of of Cis called the relative interior, and is denoted by ri(C). Definition N.7 (Ellipsoid).For a given P∈Sn ++,c∈Rn,r∈R>0we denote by E(P, c, r). ={x∈Rn:kx−ck2 P≤r2} the ellipsoid centered at cwith radius rand geometry determined by P. Definition N.8 (Admissible invariant set).Consider a discrete-time autonomous system z(t+ 1) = f(z(t)), where tis the current discrete-time instant and f:Rn→Rndescribes the system dynamics, that is subject to z(t)∈ Z ⊆ Rn,∀t. An admissible invariant set of the system is a set Ω ⊆ Z such that z(t)∈Ω =⇒z(t+ 1) ∈Ω. Functions We use the standard notation for the derivative and partial derivative of a differentiable real-valued function f. That is, df dx is the derivative of f(x) with respect to x(we may also use the shorthand f0) and ∂f ∂x is the partial derivative v
4Chapter 1. Introduction pressed sparse column/row or dictionary of keys). Instead, we only need to store the repeating submatrices, not needing to store the information of their location. The developed solvers have been made available in an open-source Matlab toolbox SPCIES [4]. The toolbox takes the model of the system and parameters of the MPC controller, and automatically generates the library-free code of the solver for its implementation in the embedded system of choice. Additionally, we present a novel MPC formulation, which we call harmonic MPC and label by HMPC, which is particularly suitable for its use with small prediction horizons, since it exhibits a larger domain of attraction that reduces the feasibility and performance issues that some MPC formulations display when the prediction horizon is not large enough. Due to its good performance with small prediction horizons, this MPC formulation may be suitable for its use in embedded systems. Part II is divided into the following chapters Chapter 4 introduces the topic of this part of the dissertation. We start by describing the problem formulation in Section 4.1 and briefly describing MPC in Section 4.2. Finally, in Section 4.3 we present a non-exhaustive review of the literature in the topic of the implementation of MPC in embedded systems. Chapter 5 presents our approach to sparsely solving MPC formulations using first order methods. We start by presenting some rather straightforward sparse algorithms for solving specific QP problems and systems of equations in Section 5.1. Section 5.2 and Section 5.3 show how we solve QP problems using two of the first order methods presented in Part I of this dissertation. This approach to solving QP problems, along with the algorithms presented in Section 5.1, will be used in Sections 5.4 to 5.6 to develop sparse solvers for the MPC formulations presented in each one of these sections. Section 5.7 presents three systems used as test benches in the numerical results presented in Section 5.8, where we compare our proposed solvers between each other and against other alternatives from the literature. Additionally, we show the results of incorporating the restart schemes of Part I to some of the proposed solvers. Chapter 6 presents the harmonic MPC formulation, starting with its description in Section 6.1. Sections 6.2 and 6.3 prove the recursive feasibility and asymptotic stability of the proposed formulation, respectively. A discussion related to the design of one of its main ingredients is presented in Section 6.4. Finally, we showcase its potential advantages in Section 6.5, where we compare its performance and domain of attraction against other MPC formulations with a particular focus on using small prediction horizons.
1.3. Publications 5 1.3 Publications The results shown in this dissertation are supported by several journal and congress publications, some of which are currently under review. The novel restart schemes and proofs presented in Part I can be found in: [5] T. Alamo, P. Krupa, and D. Limon, “Restart FISTA with global linear convergence,” in Proceedings of the European Control Conference (ECC), pp. 1969–1974, IEEE, 2019. Extended version available at arXiv:1906.09126. [6] T. Alamo, P. Krupa, and D. Limon, “Gradient based restart FISTA,” in Proceedings of the Conference on Decision and Control (CDC), pp. 3936–3941, IEEE, 2019. [7] T. Alamo, P. Krupa, and D. Limon, “Restart of accelerated first order methods with linear convergence for non-strongly convex optimization,” arXiv preprint: 2102.12387, 2021. Submitted to Transactions on Automatic Control. The sparse solvers for the implementation of MPC in embedded systems presented in Chapter 5, as well as some additional numerical results to the ones presented in this dissertation, can be found in: [8] 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. [9] P. Krupa, D. Limon, and T. Alamo, “Implementation of model predictive controllers in programmable logic controllers using IEC 61131-3 standard,” in Proceedings of the European Control Conference (ECC), pp. 288–293, IEEE, 2018. [10] 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,” arXiv preprint: 2008.09071, 2020. Submitted to Transactions on Control Systems Technology. [11] P. Krupa, J. Camara, I. Alvarado, D. Limon, and T. Alamo, “Real-time implementation of MPC for tracking in embedded systems: Application to a two-wheeled inverted pendulum,” in Proceedings of the European Control Conference (ECC), 2021. Preprint available at arXiv:2103.14571. [12] P. Krupa, R. Jaouani, D. Limon, and T. Alamo, “A sparse ADMM-based solver for linear MPC subject to terminal quadratic constraint,” arXiv preprint: 2105.08419, 2021. Submitted to Automatica.
6Chapter 1. Introduction The harmonic MPC formulation presented in Chapter 6 can be found in: [13] P. Krupa, D. Limon, and T. Alamo, “Harmonic based model predictive control for set-point tracking,” IEEE Transactions on Automatic Control, 2020. [14] P. Krupa, M. Pereira, D. Limon, and T. Alamo, “Single harmonic based model predictive control for tracking,” in Proceedings of the Conference on Decision and Control (CDC), pp. 151–156, IEEE, 2019. Other articles related to this dissertation but not directly discussed in it are: [15] P. Krupa, N. Saraf, D. Limon, and A. Bemporad, “PLC implementation of a real-time embedded MPC algorithm based on linear input/output models,” IFAC-PapersOnLine, vol. 53, no. 2, pp. 6987–6992, 2020. [16] I. Alvarado, P. Krupa, D. Limon, and T. Alamo, “Tractable robust MPC design based on nominal predictions,” arXiv preprint: 2104.06088, 2021. Submitted to Journal of Process Control. [17] P. Krupa, C. Danielson, C. Laughman, S. A. Bortoff, D. J. Burns, S. Di Cairano, and D. Limon, “Modelica implementation of centralized MPC controller for a multi-zone heat pump,” in Proceedings of the European Control Conference (ECC), 2019.
Part I Restart schemes for accelerated first order methods
9 Chapter 2 Preliminaries: Accelerated first order methods First order methods (FOM) are iterative numerical methods for solving optimization problems that only require knowledge of the gradient/subgradient (but not Hessian) of the objective function [18, 19]. The origin of these methods can be traced back to the works of Cauchy in 1847 [20]. Since then, many different FOMs have been proposed and developed [18, 19, 21, 22]. These methods have been used in a wide variety of fields, including machine learning and model predictive control [18, 23]. The rising number of fields and applications that deal with large optimization problems has increased the interest in these methods due to their lower computational cost per iteration and memory requirements when compared to other optimization methods, such as interior point or active set methods. As a gentle introduction to the topic, the most most well know FOM is the gradient descent method. Consider the optimization problem min z∈Rnzf(z),(2.1) where the convex function f:Rnz→Ris Lipschitz continuously differentiable. Then, the gradient descent method applies, for k≥0, the recursion zk+1 =zk−ρk∇f(zk) starting at an initial point z0∈Rnzand where ρk∈R>0is the step size. If ρkis chosen appropriately, the above recursion converges to an optimal solution z∗of (2.1) as k→ ∞ [18]. In 1983, Yurii Nestevor proposed a variation of the gradient descent method known as the accelerated gradient descent method [24]. This new method provided a O(1/k2) convergence rate (in terms of objective function value), which is the optimal convergence rate for the class of optimization problems (2.1); a significant improvement over the O(1/k) rate of the gradient descent method. The
10 Chapter 2. Preliminaries: Accelerated first order methods acceleration is obtained by providing a certain “momentum” to the algorithm in the form of an overrelaxation step. In particular, the iteration of the accelerated gradient descent method are, for k≥0, zk+1 =yk−ρk∇f(zk), tk+1 =1 21 + q1+4t2 k, yk+1 =zk+tk−1 tk+1 (zk−zk−1), starting at an initial point y0∈Rnzand taking t0= 1. The idea behind this method has been applied to many other FOMs, ranging from methods for smooth convex optimization [24, 25], to methods for composite non-smooth convex optimization [22, 26, 27, 28]. These methods are referred to as accelerated, of fast, first order methods (AFOM). This chapter introduces a few FOMs and AFOMs. In particular, we present two popular FOMs: the proximal gradient method and the Alternating Direction Method of Multipliers (ADMM); and the accelerated variant of the proximal gradient method, which we label as FISTA. Additionally, we present a variation of ADMM known as extended ADMM (EADMM), and a variant of FISTA known as monotone FISTA (MFISTA). We show the pseudocode of each method and some of the their properties (in particular, those that will be of use for future developments throughout this dissertation). The properties we show in this chapter are for the most part well established in the literature. However, we include our own proofs for some of them. We do this to show that the properties we cite are still applicable to our particular problem formulation and/or assumptions, which are not always the same as the ones that can be found in the references provided. We refer the reader to [18, 19, 22, 29, 30] for a few references in the field of first order methods, and to [1, 23] for general references on convex optimization. 2.1 First order methods This section describes two popular FOMs: the proximal gradient method and the alternating direction method of multipliers (ADMM). Additionally, we describe an extension of ADMM to composite convex optimization problems with three functions in the objective function. 2.1.1 Proximal gradient method The proximal gradient method is a FOM for convex composite optimization problems that is based on performing a proximal operator at each iteration of the algorithm [18, 29]. We follow, however, an equivalent description of the method, in which we make use of the so called composite gradient mapping [27, § 2], which we formerly define further ahead. For an in depth explanation of the proximal
2.1. First order methods 11 gradient method and its different variants, we refer the reader to [18, § 10.2], and, as a more accessible reference, to [29]. Consider the convex composite optimization problem f∗= min z∈Z {f(z). = Ψ(z) + h(z)},(2.2) under the following assumption. Assumption 2.1. We assume that: (i) The function Ψ ∈Γ(Rnz) may be non-smooth. (ii) h:Rnz→Ris a smooth differentiable convex function. That is, there exists R∈Snz ++ such that h(z)≤h(y) + h∇h(y), z −yi+1 2kz−yk2 R,∀(z, y)∈Rnz×Rnz.(2.3) (iii) Z ⊆ Rnzis a non-empty closed convex set. (iv) Problem (2.2) is solvable. That is, there exists z∗∈ Z ∩dom(Ψ) such that f∗=f(z∗) = inf z∈Z f(z). We denote by Ωf. ={z∈Rnz:z∈ Z, f(z) = f∗}the optimal set of problem (2.2), which, due to Assumtion 2.1.(iv) is non-empty. It is well known that this set is a singleton if fis strictly convex. Remark 2.2. It is standard to write the smoothness condition (2.3) as the well known descent lemma h(z)≤h(y) + h∇h(y), z −yi+L 2kz−yk2 S,∀(z, y)∈Rnz×Rnz,(2.4) where L∈R>0provides a bound on the Lipschitz constant of ∇h[3, § 2.1] and S∈Snz ++ is usually taken as the identity matrix. However, since L 2kz−yk2 S=1 2kz−yk2 LS,∀(z, y)∈Rnz×Rnz, we have that (2.3) implies (2.4) if we take R=LS. We use expression (2.3) because it simplifies the algebraic expressions and results we present throughout this dissertation, but analogous results can be obtained if the standard L-smoothness definition (Definition N.20) and descent lemma are used instead of (2.3). Given y∈Rnz, one could use the local information given by ∇h(y) to minimize the value of faround y. Under Assumtion 2.1, this can be done obtaining the optimal solution of the strictly convex optimization problem min z∈Z Ψ(z) + h∇h(y), z −yi+1 2kz−yk2 R.
12 Chapter 2. Preliminaries: Accelerated first order methods The solution of this optimization problem leads to the notion of the composite gradient mapping [27, § 2], which constitutes a generalization of the gradient mapping (see [3, § 2.2.3] for the case when Ψ(·) = 0 and [26] for the case when Z=Rnz). Definition 2.3 (Composite gradient mapping).Let f∈Γ(Rnz) be the composition of two functions f= Ψ + h, and let Z ⊆ Rnz. Let Assumtion 2.1 hold. We denote the following mappings Tf,Z R:Rnz→Rnzand Gf,Z R:Rnz→Rnz: Tf,Z R(y) = arg min z∈Z Ψ(z) + h∇h(y), z −yi+1 2kz−yk2 R,(2.5a) Gf,Z R(y) = Ry−T f,Z R(y).(2.5b) When the identities of f,Rand Zare clear from the context, we will often omit the superscripts and subscript and simply write T(y) and G(y). The following proposition states the uniqueness of the composite gradient mapping under Assumtion 2.1. This result is well known in the literature (see, for instance, [18, § 6.1] for an analogous result). However, we include a proof of the proposition for our particular notation and description of problem (2.2). Proposition 2.4 (On the unique solution of the composite gradient mapping). Let f∈Γ(Rnz) be the composition of two functions f= Ψ + h, and let Z ⊆ Rnz. Let Assumtion 2.1 hold. Then, Tf,Z R(y) is a singleton for any y∈Rnz. Proof: Let ΨZ= Ψ + IZ. Then, after some simple algebraic manipulations and cancellation of constant terms, we have that problem (2.5) can be equivalently written as Tf,Z R(y) = arg min z∈RnzΨZ+1 2kz−y−R−1∇h(y)k2 R.(2.6) From Assumtion 2.1.(iv), we have that dom(ΨZ) = dom(Ψ) ∩ Z 6=∅. Moreover, since both Ψ ∈Γ(Rnz) (Assumtion 2.1.(i)) and IZ∈Γ(Rnz) (since Assumtion 2.1.(iii) states that Zis non-empty), we conclude that ΨZ∈Γ(Rnz). Therefore, (2.6) is a strongly convex function that is not everywhere infinite. As such, it is well known that it has a unique minimizer for any y∈Rnz. Equation (2.6) shows that the composite gradient mapping is closely related to the proximal operator. Indeed, if R=LInzfor some L∈R>0, then by definition of the proximal operator (Definition N.24), we have that (2.6) is equivalent to proxL−1ΨZ(y−R−1∇h(y)). In the context of optimal gradient methods, it is assumed that the computation of T(y) is cheap. This is the case when Zis a simple set (a box, Rnz, etc.), R∈Dnz ++, and Ψ(·) is a separable function. For example, in the well known Lasso optimization problem, the computation of T(y) resorts to the computation
2.1. First order methods 13 of the shrinkage operator [26]. See [18, § 6], [29, § 6],[31], or [32, § 28] for numerous examples in which the computation of the composite gradient mapping is simple. The following proposition gathers some well-known properties of the composite gradient mapping and the dual norm kG(·)kR−1. We include the proof for completeness. Proposition 2.5. Consider problem (2.2) and let Assumtion 2.1 hold. Then: (i) For every y∈Rnand z∈ Z, f(T(y)) −f(z)≤ hG(y),T(y)−zi+1 2kG(y)k2 R−1(2.7a) =hG(y), y −zi− 1 2kG(y)k2 R−1(2.7b) =−1 2kT(y)−zk2 R+1 2ky−zk2 R.(2.7c) (ii) For every y∈ Z,1 2kG(y)k2 R−1≤f(y)−f(T(y)) ≤f(y)−f∗. Proof: Let hy:Rnz→Rbe the given by hy(z) = h∇h(y), z −yi+1 2kz−yk2 R, ΨZ:Rnz→Rbe given by ΨZ(z) = Ψ(z) + IZ(z), and Fy:Rnz→Rbe given by Fy(z)=ΨZ(z) + hy(z). From Proposition 2.4, we have that, for any given y∈Rnz,T(y) is the unique minimizer of problem (2.5a). Therefore, Ψ(T(y)) + hy(T(y)) ≤Ψ(z) + hy(z),∀z∈ Z, which, since T(y)∈ Z, implies ΨZ(T(y)) + hy(T(y)) ≤ΨZ(z) + hy(z),∀z∈Rnz. From this inequality we have that Fy(T(y)) ≤Fy(z),∀z∈Rnz, which due to the definition of the subdifferential (Definition N.18) implies 0nz∈∂Fy(T(y)).(2.8) From Assumtion 2.1.(iv), we have that dom(ΨZ) = dom(Ψ) ∩Z 6=∅. Moreover, since both Ψ ∈Γ(Rnz) (due to Assumtion 2.1.(i)) and IZ∈Γ(Rnz) (since Assumtion 2.1.(iii) states that Zis non-empty), we conclude that ΨZ∈Γ(Rnz). Additionally, since hyis a continuous real-valued function in Rnz, we have that it is closed (see [1, Proposition 1.1.3]) and that its domain is Rnz. Therefore, Fyis the sum of two closed convex functions satisfying ri(dom(ΨZ)) ∩ri(dom(hy)) 6=∅, which along [1, Proposition 5.4.6], lets us conclude that ∂Fy(T(y)) = ∂ΨZ(T(y)) + ∂hy(T(y)), where the subdifferential ∂hyof the differentiable function hyis given by ∂hy(T(y) = ∇hy(T(y)) = ∇h(y) + R(T(y)−y).
20 Chapter 2. Preliminaries: Accelerated first order methods Proof: Claim (i) is trivial, since tkis defined as the roots of t2 k−tk−t2 k−1= 0 for all k≥1. Claim (ii) is proven by induction on k. The claim is trivially satisfied for k= 0. Suppose that the claim is satisfied for k−1, i.e. tk−1≥k+2 2. Then, tk=1 21 + q1+4t2 k−1≥1 21 + q4t2 k−1=1 2+tk−1≥1 2+k+ 1 2=k+ 2 2, which shows that it is also satisfied for k. Proposition 2.12. Consider problem (2.2) and let Assumtion 2.1 hold. Let {zk}≥0and {yk}≥0be the sequences generated by FISTA (Algorithm 4) for problem (2.2) starting at z0=y0=Tf,Z R(z), where z∈Rnzis given. Then: (i) f(zk)−f∗≤2 (k+ 1)2kz0−πR Ωf(z0)k2 R,∀k≥1. (ii) kGf,Z R(yk)kR−1≤4 k+ 2kz0−πR Ωf(z0)kR,∀k≥0. Proof: For the sake of simplicity, we use the shorthand notations ¯z0. =πR Ωf(z0), Gk. =Gf,Z R(yk), y+ k . =Tf,Z R(yk), δfk. =f(zk)−f∗,δzk. =zk−¯z0and δyk. =yk−¯z0, where the subscripts may change, e.g., y+ k+1 . =Tf,Z R(yk+1). We start by proving claim (i). Particularizing inequality (2.7c) from Proposition 2.5.(i) to y0∈ Z ⊆ Rnzand ¯z0∈Ωf⊆ Z, we obtain f(y+ 0)−f(¯z0)≤ −1 2ky+ 0−¯z0k2 R+1 2ky0−¯z0k2 R. From steps 1 and 4 of the algorithm we have z0=y0and z1=y+ 0. Furthermore, by definition of ¯z0we have f(¯z0) = f∗. Therefore we can rewrite the previous inequality as f(z1)−f∗≤ −1 2kz1−¯z0k2 R+1 2kz0−¯z0k2 R≤1 2kz0−¯z0k2 R,(2.15) which proves claim (i) for k= 1. We now proceed to prove the claim for k≥2. From step 4 of Algorithm 4 we have zk=y+ k−1,∀k≥1.(2.16) Therefore, inequality (2.7b) from Proposition 2.5.(i) leads to f(z)≥f(zk+1) + 1 2kGkk2 R−1−hGk, yk−zi,∀z∈ Z,∀k≥1. We notice that, by construction, zk∈ Z,∀k≥1. Thus, particularizing the previous inequality to z=zkand to z= ¯z0∈Ωf⊆ Z, we obtain the inequalities: f(zk)≥f(zk+1) + 1 2kGkk2 R−1−hGk, yk−zki,∀k≥1, f(¯z0)≥f(zk+1) + 1 2kGkk2 R−1−hGk, yk−¯z0i,∀k≥1,
2.2. Accelerated first order methods 21 which can be rewritten as δfk−δfk+1 ≥1 2kGkk2 R−1−hGk, δyk−δzki,∀k≥1,(2.17a) −δfk+1 ≥1 2kGkk2 R−1−hGk, δyki,∀k≥1.(2.17b) We introduce now the auxiliary nomenclature Γkdefined as Γk. =t2 k−1δfk−t2 kδfk+1,∀k≥1, which, making use of Lemma 2.11.(i) allows us to write Γk= (t2 k−tk)δfk−t2 kδfk+1 = (t2 k−tk)(δfk−δfk+1)−tkδfk+1,∀k≥1.(2.18) In view of Lemma 2.11.(ii),tk≥1,∀k≥0. This implies that we can replace the terms δfk−δfk+1 and −δfk+1 in inequality (2.18) by the lower bounds given in inequalities (2.17). Doing so, we obtain: Γk≥(t2 k−tk)1 2kGkk2 R−1−hGk, δyk−δzki+tk1 2kGkk2 R−1−hGk, δyki =t2 k 2kGkk2 R−1−hGk, t2 k(δyk−δzk) + tkδzki,∀k≥1.(2.19) Step 6 of the algorithms can be rewritten as δyk−δzk=tk−1−1 tk (δzk−δzk−1),∀k≥1, which, introducing the notation sk. =δzk−1+tk−1(δzk−δzk−1), leads to sk−δzk=δzk−1+tk−1(δzk−δzk−1)−δzk= (tk−1−1)(δzk−δzk−1) =tk(δyk−δzk),∀k≥1.(2.20) Using the definition of skalong with (2.20), we show that Gkcan be written in terms of skand sk+1: tkGk (∗) =tkR(yk−zk+1) = tkR(δyk−δzk+1) =tkR(δyk−δzk+δzk−δzk+1) =R(sk−δzk+tk(δzk−δzk+1)) =R(sk−sk+1),∀k≥1,(2.21) where (∗) is due to (2.16). We can now keep developing (2.19) as follows Γk≥1 2ktkGkk2 R−1−hGk, tk(sk−δzk) + tkδzki=1 2ktkGkk2 R−1−htkGk, ski (2.21) ≥1 2kR(sk−sk+1)k2 R−1−hR(sk−sk+1), ski =1 2ksk+1 −skk2 R+hR(sk+1 −sk), ski =1 2k(sk+1 −sk) + skk2 R−1 2kskk2 R=1 2ksk+1k2 R−1 2kskk2 R,∀k≥1.
22 Chapter 2. Preliminaries: Accelerated first order methods Thus, for every k≥1, Γk=t2 k−1δfk−t2 kδfk+1 ≥1 2ksk+1k2 R−1 2kskk2 R. Equivalently, t2 kδfk+1 +1 2ksk+1k2 R≤t2 k−1δfk+1 2kskk2 R,∀k≥1. Since this inequality holds ∀k≥1, we can apply it in a recursive way to obtain t2 kδfk+1 +1 2ksk+1k2 R≤t2 0δf1+1 2ks1k2 R =δf1+1 2kδz0+t0(δz1−δz0)k2 R =δf1+1 2kz1−¯z0k2 R (2.15) ≤1 2kz0−¯z0k2 R,∀k≥1. Therefore, t2 k(f(zk+1)−f∗) + 1 2ksk+1k2 R≤1 2kz0−¯z0k2 R,∀k≥1, from where, making use of Lemma 2.11.(ii), we finally conclude that f(zk+1)−f∗≤kz0−¯z0k2 R 2t2 k≤2kz0−¯z0k2 R (k+ 2)2,∀k≥1. Next, we prove claim (ii). We start by proving the claim for k= 0. From (2.15) we derive kz1−¯z0kR≤ kz0−¯z0kR.(2.22) Thus kG(y0)kR−1=kR(y0−y+ 0)kR−1=ky0−y+ 0kR=kz0−z1kR =kz0−¯z0+ ¯z0−z1kR≤ kz0−¯z0kR+kz1−¯z0kR ≤2kz0−¯z0kR. We now prove the claim for k > 0. From s1=δz0+t0(δz1−δz0) = z1−¯z0, along with (2.22), we derive ks1kR=kz1−¯z0kR≤ kz0−¯z0kR, which in addition to the following inequality that can be deduced from (2.15): ksk+1kR≤ kz0−¯z0kR,∀k≥1, leads to kskkR≤ kz0−¯z0kR,∀k≥1. From here we derive, for every k≥1, ksk+1−skkR≤ ksk+1kR+kskkR≤ kz0−¯z0kR+kz0−¯z0kR= 2kz0−¯z0kR.(2.23)
2.2. Accelerated first order methods 23 Algorithm 5: MFISTA Require: z∈dom f,∈R>0 1y0←z,z0←z,t0←1, k←0 2repeat 3k←k+ 1 4vk← Tf,Z R(yk−1) 5tk←1 21 + q1+4t2 k−1 6zk← vkif f(vk)≤f(zk−1) zk−1otherwise 7yk←zk+tk−1 tk (vk−zk) + tk−1−1 tk (zk−zk−1) 8until kG(zk)kR−1≤ Output: ˜z∗←zk From (2.21) we have Gk=1 tk R(sk−sk+1),∀k≥1. Therefore, taking into account Lemma 2.11.(ii), kGkkR−1=1 tkksk−sk+1kR (2.23) ≤2 tkkz0−¯z0kR≤4 k+ 2kz0−¯z0kR,∀k≥1. 2.2.2 A monotone variant of FISTA This section presents a monotone version of FISTA, which we label by MFISTA, and which is presented in [42, § V.A]. FISTA is not a monotone algorithm, in that the function values {f(zk)}≥0generated by Algorithm 4 are not necessarily nonincreasing. However, a small adjustment can be made to it to provide a monotone behaviour. Algorithm 5 shows the MFISTA algorithm applied to problem (2.2) under Assumption 2.1 for an initial point z∈dom(f) and an exit tolerance ∈R>0. It returns a suboptimal solution ˜z∗of problem (2.2), where the suboptimality is determined by . We note that Remark 2.10 also applies to Algorithm 5. Note that MFISTA and FISTA both require a single computation of the composite gradient mapping Tf,Z R(yk−1) and a few vector-vector operations (with MFISTA requiring a few more than FISTA). The main difference between the two, computationally speaking, is that MFISTA requires the two function evaluations f(vk) and f(zk−1) at each iteration, while FISTA requires none. The interesting aspect of MFISTA is that it shares FISTA’s convergence rate result given in Proposition 2.12.(i), as stated in the following proposition.
24 Chapter 2. Preliminaries: Accelerated first order methods Proposition 2.13. Consider problem (2.2) and let Assumtion 2.1 hold. Let {zk}≥0be the sequence generated by MFISTA (Algorithm 5) for problem (2.2) starting at z∈dom(f). Then, f(zk)−f∗≤2 (k+ 1)2kz0−πR Ωf(z0)k2 R,∀k≥1. Proof: The proof of this proposition requires modifications to the proof given for FISTA (see Proposition 2.12). We refer the reader to [42] for its proof. 2.3 Conclusions In this chapter we have presented some well known FOMs and AFOMs. As we discussed previously, the popularity of these methods has increased in recent years do to the rising number of applications that deal with large optimization problems. In this case, even though FOMs typically require more iterations to solve than methods that make use of the Hessian of the optimization problem, their lower computational cost per iteration may make them a more suitable choice. However, these methods are not without issue. In addition to the fact that they may require a large number of iterations to converge, there are other aspects that have drawn the attention of the scientific community. One of the main challenges is the study and improvement of their theoretical and practical convergence properties. As an example, numerous results address the problem of selecting the penalty parameter of ADMM [22, 43, 44]. In the following chapter we will address one of these challenges: improvements of the convergence of AFOMs related to the undesirable oscillating behaviour that they display when close to the optimal solution.
25 Chapter 3 Restart schemes for accelerated first order methods One of the drawbacks of accelerated first order methods is that they may suffer from undesirable oscillatory behavior which slows down their convergence [45]. The intuitive reason behind this is that AFMOs can be thought of as momentum driven, in the sense that each iterations depends on the previous ones with an added momentum that is increased over time. This momentum, coupled with the fact that most AFOMs are non-monotone in the objective function value, leads to the appearance, in many applications, of a periodic oscillation when close to the optimal solution. To illustrate this issue, let us take a look at a simple example. Example 3.1. Consider the quadratic programming (QP) problem min z∈R2{f(z). =1 2z>Hz +q>z},(3.1) where H= diag(0.5,1) and q=−(0.1,1). It is well known that the solution of this problem is given by z∗=−H−1q= (0.2,1). We solve problem (3.1) using FISTA algorithm (Algorithm 4) with the initial condition z=−(2,5), and taking R= 100I2and = 10−6. Remark 3.2. We note that our selection of Rin Example 3.1, while it does satisfy Assumtion 2.1.(ii) for the given smalll-sized optimization problem, could have easily been selected to provide a much better performance of FISTA. For larger optimization problems, however, it is common for Rto be overestimated, since the computation of a better value is not worth the computational burden. For the purposes of this example, it was selected this way to enhance the oscillatory behavior of the algorithm, thus providing visually obvious oscillations. We remark that the oscillatory behavior would be present, in a greater or lesser extent, when solving most optimization problems if the exit tolerance is sufficiently small. In Section 3.3 we show examples with larger-sized optimization problems.
26 Chapter 3. Restart schemes for accelerated first order methods -0.5 0 0.5 1 0 0.5 1 1.5 2 (a) Trajectory of zk.(b) Convergence to optimal objective function value. Figure 3.1: Example 3.1. FISTA without a restart scheme. Figure 3.1 shows the result of solving Example 3.1 using FISTA (Algorithm 4). Figure 3.1a shows the (incomplete) trajectory of the iterates zkproduced by FISTA, whereas Figure 3.1b shows the difference f(zk)−f∗in a logarithmic scale. Notice that during the initial iterations f(zk) is monotonically converging towards the optimal solution, as can be seen in Figure 3.1b. However, it soon becomes non-monotone with the appearance of the oscillations. The general trend still converges towards the optimal solution, but at a slower rate. The algorithm converged after 853 iterations to the given tolerance. This issue can be addressed by means of a restart scheme. Intuitively, a restart scheme is a procedure in which an iterative method is stopped when a certain condition, which we label the restart condition, is met, and is then restarted using the last value of the iterates as the new initial condition. This process continues until a terminal exit condition is satisfied. In the context of AFOMs, and as a more formal description of a restart scheme, we consider the following. Let Abe a non-specific AFOM that starts at a given initial point z0, and generates a sequence {zk}≥0that converges towards the optimal solution z∗of the optimization problem as k→ ∞. The algorithm terminates at a finite iterate kout when a certain condition Ecis met. In a non-restarted paradigm, this restart condition will be a certificate for the degree of suboptimality of zout . =zkout . For instance, the exit condition of FISTA (Algorithm 4) is given by Ec(zk) = (true if kG(zk)kR−1≤, false otherwise, where ∈R>0is the given exit tolerance. The exit condition Ecis therefore a map to {true,false}in relation to the current state of the algorithm. If the exit condition Ecis evaluated to true, then the algorithm terminates, returning zout ←zkand kout ←k. We will write this as [zout, kout]← A(z, Ec). A restarted paradigm follows the procedure shown in Algorithm 6. The procedure uses an AFOM algorithm A; an initial condition r0∈ C ⊆ Rnz, where Cwill depend on the algorithm A; a restart condition Ec; and a terminal exit
3.1. A brief review of the literature 27 Algorithm 6: General restart procedure for AFOMs. Require: r0∈ C,Ec,Et 1k←0, j←0 2repeat 3j←j+ 1 4[rj, kj]← A(rj−1, Ec) 5k←k+kj 6Evaluate Et; 7until Et=true Output: rout ←rj,jout ←j,kout ←k condition Et. The algorithm makes successive calls to A, generating a sequence {rj}≥0of restart points, where jis the restart counter. At each iteration j(with the exception of the first iteration) algorithm Ais given the value it returned in its previous call and is executed until the restart condition Ecis satisfied. The terminal exit condition Etshould be a measure of the optimality of rj, whereas the restart condition Eccan be any condition; although, obviously, it will be chosen to improve the (practical) convergence of this scheme when compared with the non-restarted variant. The algorithm returns, in rout, the last rjprovided by A, the total number of restarts jout, and the total number kout of iterations performed by Aamong all its calls. In the following, we will use the letters rand jto refer to the restart points and restart counter, respectively, and the letters z and kto refer to the inner iterates of A. 3.1 A brief review of the literature Various restart schemes for AFOMs have been presented in the literature [27, 45, 46, 47, 48, 49, 50, 51], each considering a particular class of optimization problem and with different assumptions on the knowledge (or lack thereof) of the parameters that characterize it. Additionally, each one derives convergence properties of their restart scheme under certain assumptions. These convergence properties are often based on the notion of the iteration complexity of the algorithm, which is the upper bound of the number of iterations required to find an -accurate optimal solution. Typically, and unless we explicitly state otherwise, this is expressed as the number of iterations krequired to reach a suboptimal solution f(zk) satisfying f(zk)−f∗≤, but other measures of suboptimality may be used. Complexity bounds have been derived for many first order methods and their accelerated variants. For instance, when applied to an L-smooth convex optimization problem, the projected gradient method attains a complexity O(L/), whereas
28 Chapter 3. Restart schemes for accelerated first order methods the fast projected gradient method has the much better complexity O(pL/) [3, 26, 52]. If, in addition, the problem is µ-strongly convex, then the iteration complexities are O(L/µ log(1/)) for the projected gradient method and O sL µlog 1 !(3.2) for the fast projected gradient method [3, 27, 49]. The log(1/) term is particularly relevant here, since it is indicative of a linear convergence rate, i.e., when the iterates satisfy an expression of the form f(zk)−f∗≤Cρk, for some C∈Z>0and ρ∈(0,1) [53, § 3.1]. Furthermore, the iteration complexity (3.2) is the optimal one for this class of optimization problems using a first order black box oracle [3], meaning that a better iteration complexity cannot be attained (up to a constant, that is) with a first order method that does not take advantage of particularities of the problem structure. One of the main questions when dealing with restart schemes is whether or not iteration complexities can be derived; in particular, which of the above and under what assumptions. Since restart schemes will be applied to AFOMs, we are particularly interested in determining under what conditions (if any), does a restart scheme attain the optimal iteration complexity (3.2). 3.1.1 Fixed-rate restart schemes The most simple class of restart schemes are what we refer to as fixed-rate restart schemes, in which the AFOM is restarted every time it performs a predetermined fixed number of iterations. That is, the restart condition is given by Ekm(k) = (true if k≥km false otherwise, where km∈Z>0is fixed. In [27, § 5.1], a fixed-rate restart scheme is presented which considers the class of optimization problems minz∈Z f(z), where Z ⊆ Rnzis a closed convex set and f:Rnz→Ris a µ-strongly convex and L-smooth function. By choosing the restart rate appropriately, this restart scheme recovers the optimal iteration complexity (3.2) for this class of optimization problem. Similar (and oftentimes equivalent) fixed-rate restart schemes are presented in other publications, such as in [48, § 3.2], [45, § 3.1]. Another particularly relevant fixed-rate restart scheme was presented in [46, § 5.2.2]. The interest of this scheme is that it proves linear convergence of the restarted scheme for non-necessarily strongly convex optimization problems. In particular, the paper considers the class of optimization problems f∗= min z∈Z f(z),(3.3)
3.1. A brief review of the literature 29 where Z ⊆ Rnzis a closed convex set and f:Rnz→Ris a L-smooth convex function satisfying the following condition, known as the quadratic functional growth condition (see [46, Definition 4]), for some µ > 0. Definition 3.3 (Quadratic functional growth).Consider problem (3.3), where Z ⊆ Rnzis a closed convex set and f:Rnz→R, and assume that its optimal set Ωfis nonempty and that f∗is finite. We say that fhas a quadratic functional growth on Zif there exists a constant µ∈R>0such that f(z)−f∗≥µ 2kz−πΩf(z)k2 2, z ∈ Z.(3.4) Remark 3.4. Inequality (3.4) is satisfied, at least locally, for a large class of not necessarily strongly convex nor strictly convex functions [46, 54]. We note that it is always satisfied if fis strongly convex. The scheme presented in [46, § 5.2.2] can be implemented in two different ways. The first is by performing the fixed-rate restart given by Ek∗(k) = (true if k≥k∗, false otherwise, where k∗=&2esL µ'. The second alternative is to restart Awhenever the following restart condition E∗ f:Rnz→ {true,false}is satisfied: E∗ f(zk) = true if f(zk)−f∗≤f(z0)−f∗ e2, false otherwise. (3.5) We note that this second approach, while equivalent to the first in the sense that the same convergence results are obtained, is not (strictly speaking) a fixed-rate restart scheme. There are two main problems with fixed-rate restart schemes. First, they require knowledge of parameters that are not typically known beforehand or that are often computationally expensive to compute, such as the parameters that characterize the smoothness, the strong convexity, the quadratic functional growth or the optimum value f∗. Second, the restart rate is derived from global parameters of the optimization problem, and may thus result in an inappropriate restart rate in better conditioned regions [45, § 3.1]. Example 3.5. As an example of the use of a fixed-rate restart scheme, Figure 3.2 shows the use of the restart scheme from [46, § 5.2.2] using the exit condition (3.5) applied to FISTA to solve Example 3.1, along with the result of the non-restarted
36 Chapter 3. Restart schemes for accelerated first order methods Algorithm 7: Restart FISTA based on objective function values Require: r0∈ Z, > 0 1n0←0, j←1k←0 2[r1, n1]←FISTA r0, Ee f(nj−1) 3repeat 4j←j+ 1 5[rj, nj]←FISTA rj−1, nj−1, Ee f(nj−1) 6k←k+nj 7if f(rj−1)−f(rj)>f(rj−2)−f(rj−1) ethen 8nj←2nj−1 9end if 10 until kG(rj)kR−1≤ Output: rout ←rj,jout ←j,kout ←k Proposition 3.11.(iii) tells us that there exists a finite number of iterations after which the distance to the optimal value f∗is smaller than a fraction of the distance covered from z0to the current iterate. Furthermore, notice that j2√e+1 √µk>j2 √µk. Therefore, if k≥j2√e+1 √µk, not only is Proposition 3.11.(iii) satisfied, but also, from Proposition 3.11.(ii), we have that f(z0)−f(zk)≥0. Therefore, a fixedrate restart scheme with the rate set to km=j2√e+1 √µk>j2 √µkwould provide linear convergence. However, this restart scheme would suffer from the same issues of other fixed-rate restart schemes: the value of µwould be required for its implementation, and the use of the global µmight lead to a worse practical convergence than if local information was used instead (see Section 3.1.1). Remark 3.12. We note that claims (ii) and (iii) of Proposition 3.11 may be satisfied for values of ksmaller that the given bounds. That is, the proposition states that the satisfaction of both inequalities is guaranteed if kis larger or equal to the given bounds, but they may also be satisfied for smaller values of k. We propose the restart scheme shown in Algorithm 7, which uses the restart condition Ee f:Rnz×Rnz×Z>0→ {true,false}given by Ee f(zk, zm, k;n) = true ⇐⇒ f(zm)−f(zk)≤f(z0)−f(zm) e(3.12a) f(zk)≤f(z0) (3.12b) k≥n, (3.12c) where m=k 2+1. For convenience, we drop the (zk, zm, k) notation, and instead simply write Ee f(n). Starting from an initial point r0∈ Z, and for a given exit tolerance ∈R>0, Algorithm 7 makes successive calls to FISTA using the above
3.2. Implementable restart schemes with linear convergence 37 exit condition. The algorithm returns an -accurate solution rout of problem (3.9), in terms of kG(rout)kR−1≤(which is a valid condition for suboptimality of rout, as shown in Proposition 2.6); the total number of restart iterations jout; and the total number of iterations kout performed by FISTA. The idea behind this restart scheme is to use an approximation of claims (ii) and (iii) of Proposition 3.11, where for claim (iii) zkis used as our best estimation of z∗and zmis the middle point of the sequence of iterates generated by the current call to FISTA. Therefore, the restart condition is only true if (i) the current iterate zkof FISTA is no worse, in terms of its objective function value, than the current initial condition z0, and (ii) if the amount that the objective function value has decreased from z0to zmis etimes smaller than the amount decreased from zmto zk. That is, (3.12a) serves as a way of detecting a degradation of the performance of FISTA, whereupon the condition is satisfied if the second half of the iterates are providing significantly less benefits that the first half, measured in terms of the decrease obtained in the objective function value; whereas (3.12b) is guaranteeing that, at the very least, the new restart point is no worse than the previous one. Remark 3.13. In fact, the imposition of (3.12b) guarantees that f(rj−1)> f(rj) provided that rj−16=z∗. This follows from noting that at every restart iteration jof Algorithm 7, step 1 of Algorithm 4 performs the assignment z0← T(rj−1). Thus, in view of Proposition 2.5.(ii) and the use of the exit condition (3.12b), we have that f(rj−1)> f(z0)≥f(zout) = f(rj)if rj−16=z∗. Therefore, in general, a net improvement is attained at every iteration j. The problem is that, due to the non-monotone nature of FISTA, the amounts f(zm)−f(zk) and f(z0)−f(zm) are not necessarily positive. However, we know from Proposition 3.11 that this will not the case if kis sufficiently large. For this reason, the restart condition also includes a minimum number of iterations nthat is adapted between calls to FISTA. In particular, at each iterate jof Algorithm 7, we set the minimum number of iterations of FISTA as nj, which is generally taken as the number of iterations performed by FISTA in iteration j−1. However, this minimum number of iterations is doubled if the condition f(rj−1)−f(rj)>f(rj−2)−f(rj−1) e(3.13) is satisfied, which, once again, is an approximation of Proposition 3.11.(iii) but this time between restart points. The idea is that njwill eventually converge to a number of iterations in which claims (ii) and (iii) of Proposition 3.11 are being satisfied by the outputs of FISTA obtained at each iteration j. If these conditions are being satisfied for the current value of nj, then it is not increased. However, if we detect that they are not satisfied, either because FISTA exits with a number of iterations njlarger than nj−1or because condition (3.13) is satisfied, then nis increased. The use of njcan be viewed as a way of finding the optimal restart rate, in the sense of the rate used by the fixed-rate restart schemes (see
38 Chapter 3. Restart schemes for accelerated first order methods -0.5 0 0.5 1 0 0.5 1 1.5 2 (a) Trajectory of zk.(b) Convergence to optimal objective function value. Figure 3.5: In red, Algorithm 7 applied to Example 3.1. The non-restarted FISTA is depicted in blue. The red circles in (a) depict each of the iterates zk. The black crosses depict each of the restart points rj. Section 3.1.1). However, we do not use global information. Instead, if the iterates a providing a good performance, then we do not increase njartificially. The main drawbacks of this restart scheme are that (i) it requires evaluating the objective function value at every iteration of FISTA, which is not necessary in the non-restarted FISTA algorithm, and (ii) the evaluation of (3.12a) requires storing the historic of the past objective function values for the current call to FISTA. This latter issue is not a significant problem in practical terms, since the use of a maximum desired number of iterations in a real implementation of iterative algorithms is commonplace. Thus, memory can be allocated for storing the past values of the objective function values even if the value of µis unknown. Example 3.14. Figure 3.5 shows the use of Algorithm 7 to solve Example 3.1 along with the result of the non-restarted FISTA algorithm. The total number of FISTA iterations is kout = 237 and the number of restarts is jout = 8. The main properties of Algorithm 7 are summarized in the following theorem. Claim (i) simply provides an upper bound of kG(rj−1)kR−1in terms of the difference f(rj−1)−f(rj), which due to (3.12b) is non-negative; claim (ii) provides an upper bound to the value of nj; and claim (iii) provides the iteration complexity of the algorithm to reach an -accurate solution of problem (3.9) in terms of kG(rout)kR−1≤. Note that the iteration complexity is given by O1 √µlog 1 ,(3.14) which we note is optimal in the sense discussed in Section 3.1 and shown in (3.2). Remark 3.15. The reader may have noted that the Lipschitz constant Ldoes not seem to be present in the iteration complexity (3.14), as is the case of the iteration complexities provided in Section 3.1 for the AFOMs. However, this is not the case, since parameter Reffectively accounts for L(see Remark 2.2). If we were to rewrite the characterization of the smoothness of husing the L-smoothness
3.2. Implementable restart schemes with linear convergence 39 condition given in Definition N.20 instead of using (2.3), then √Lwould appear in (3.14). Theorem 3.16 (Convergence properties of Algorithm 7).Consider problem (3.9) under Assumtions 2.1 and 3.9. Then, the sequences {rj}and {nj}generated by Algorithm 7 and its output kout satisfy: (i) 1 2kGf,Z R(rj−1)k2 R−1≤f(rj−1)−f(rj),∀j≥1. (ii) nj≤4√e+ 1 √µ,∀j≥0. (iii) The total number of iterations of FISTA required to attain the exit condition kGf,Z R(rout)kR−1≤is upper bounded by kout ≤16 √µlog 1 + 2(f(r0)−f∗) 2. Proof: By construction, rj−1∈ Z, for all j≥1. Therefore, we have from Proposition 2.5.(ii), that 1 2kG(rj−1)k2 R−1≤f(rj−1)−f(T(rj−1)),∀j≥1.(3.15) We also notice that rjis computed invoking FISTA algorithm using rj−1as initial condition (z=rj−1). That is, [rj, nj]←F IST A(rj−1, Ee f(nj−1)).Since the output value f(rj) is forced to be no larger than the one corresponding to z0=T(z) = T(rj−1), we have f(rj)≤f(T(rj−1)). Therefore, we obtain from inequality (3.15) that 1 2kG(rj−1)k2 R−1≤f(rj−1)−f(T(rj−1)) ≤f(rj−1)−f(rj),∀j≥1, which proves claim (i). Let us now prove that if nj−1≤4√e+1 √µ, then the value njobtained from [rj, nj]←FISTA(rj−1, Ee f(nj−1)),also satisfies nj≤4√e+ 1 √µ.(3.16) Let us denote ¯m. =j2√e+1 √µk. Since ¯m≥j2√e+1 √µk, we infer, from Proposition 3.11.(iii) that f(z¯m)−f∗≤f(z0)−f(z¯m) e. From this inequality, we obtain f(z¯m)−f(zk)≤f(z¯m)−f∗≤f(z0)−f(z¯m) e.
40 Chapter 3. Restart schemes for accelerated first order methods Therefore, the restart condition (3.12a) is satisfied for m= ¯m. Since m=k 2+1 we have m≥k 2. This means that for m= ¯m, the corresponding value for kis no larger than k≤2 ¯m= 2 2√e+ 1 √µ≤4(√e+ 1) √µ. We also notice that, in view of Proposition 3.11.(ii), the restart condition (3.12b) is satisfied for every k≥j2 √µk.Therefore, nj−1≤4√e+1 √µimplies that nj, obtained from [rj, nj]←FISTA(rj−1, Ee f(nj−1)), also satisfies (3.16). Making use of the previous discussion, we now prove claim (ii) by reduction to the absurd. Suppose that nj>4√e+ 1 √µ.(3.17) Because of the previous discussion, an taking into account that n0= 0, inequality (3.17) could only be attained by the doubling step nj= 2nj−1of the algorithm (see step 8). That is, inequality 3.17 is possible only if there is s∈Z>0>1 such that ns−1>2√e+1 √µand f(rs−1)−f(rs)>f(rs−2)−f(rs−1) e. Since [rs−1, ns−1]←FISTA(rs−2, Ee f(ns−2)),we have that rs−1is obtained from rs−2applying ns−1>2√e+1 √µiterations of FISTA algorithm. However, we have from Proposition 3.11.(iii) that this number of iterations implies f(rs−1)−f(rs)≤f(rs−1)−f∗≤f(T(rs−2)) −f(rs−1) e. From Proposition 2.5.(ii) we also have f(T(rs−2)) ≤f(rs−2), which when combined with the above inequality leads to f(rs−1)−f(rs)≤f(rs−2)−f(rs−1) e. That is, there is no doubling step if ns−1≥2√e+1 √µ, which proves claim (ii). We now show that there is a doubling step (i.e., step 8 is executed) at least every T. =log 1 + 2(f(r0)−f∗) 2 iterations of the algorithm. Suppose that there is no doubling step from iteration j=s+ 1 to j=s+T, where s∈Z>0≥1. That is, f(rj−1)−f(rj)≤f(rj−2)−f(rj−1) e,∀j∈[s+ 1, s +T].
3.2. Implementable restart schemes with linear convergence 41 From this, and claim (i) of the theorem, we obtain the following sequence of inequalities: 1 2kG(rs+T−1)k2 R−1≤f(rs+T−1)−f(rs+T)≤f(rs+T−2)−f(rs+T−1) e ≤1 eT (f(rs−1)−f(rs)) ≤1 eT (f(rs−1)−f∗) ≤1 eT (f(r0)−f∗) = 1 ellog1+ 2(f(r0)−f∗) 2m (f(r0)−f∗) ≤1 elog1+ 2(f(r0)−f∗) 2(f(r0)−f∗) = 1 1 + 2(f(r0)−f∗) 2 (f(r0)−f∗)≤2 2. We conclude that Tconsecutive iterations without doubling step implies that the exit condition kG(rs+T−1)kR−1≤is satisfied. Therefore, there must be at least one doubling step every Titerations. This implies that there exist j∈[s+1, s+T] such that f(rj−1)−f(rj)>f(rj−2)−f(rj−1) e, which, in view of step 7, implies that nj= 2nj−1. Moreover, since {nj}is a non-decreasing sequence, we get ns+T≥nj= 2nj−1≥2ns,∀s≥1. That is, ns≤ns+T 2,∀s≥1.(3.18) Let us rewrite jas j=m+nT, where 0 ≤m<T and n≥0. From the non decreasing nature of {nj}, we have that j X i=0 ni= m+nT X i=0 ni= m X i=0 ni+ n−1 X `=0 T X i=1 nm+i+`T ≤Tnm+T n X `=1 nm+`T =T n X `=0 nm+`T =T n X `=0 nj−`T .(3.19) Also, from inequality (3.18), we have nj−T≤nj 2. Using this inequality in a recursive manner we obtain nj−`T ≤1 2` nj, ` = 0, . . . , n, which along with (3.19) leads to j X i=0 ni≤T n X `=0 1 2` nj≤T∞ X `=0 1 2` nj= 2Tnj.
42 Chapter 3. Restart schemes for accelerated first order methods Algorithm 8: Gradient Based Restart FISTA Require: r0∈Rn,∈R>0 1y0← T(r0), z0← T(r0), t0←1, k←0j←0, ρ0← kG(r0)kR−1 2repeat 3k←k+ 1 4zk← T(yk−1) 5tk←1 21 + q1+4t2 k−1 6yk←zk+tk−1−1 tk (zk−zk−1) 7if kG(yk)kR−1≤ρj ethen 8j←j+ 1 9rj←yk 10 ρj← kG(rj)kR−1 11 yk← T(rj), zk← T(rj), tk←1 12 end if 13 until ρj≤ Output: zout ←zk,jout ←j,kout ←k By now making the summation up to j=jout, claim (iii) now directly follows by noting that claim (ii) states that nj≤4√e+1 √µ,∀j≥0, and by noting that kout =Pjout j=0 njif no doubling steps (step 8 of the algorithm) are performed and kout ≤Pjout j=0 njotherwise. 3.2.2 A gradient based restart scheme for FISTA This section presents a computationally cheap and simple to implement restart scheme for FISTA (Algorithm 4) that was originally presented in [6] and that exhibits linear convergence when applied to the class of optimization problems (3.9) under Assumtions 2.1 and 3.9. The scheme is based on monitoring the evolution of the composite gradient mapping Gduring the iterates of FISTA through the use of the exit condition Ee g:Rnz→ {true,false}given by Ee g(yk) = true if kG(yk)kR−1≤1 ekG(z)kR−1 false otherwise, (3.20) where ykis the iterate of FISTA computed at step 6 of Algorithm 4. That is, FISTA is restarted every time the value of the composite gradient mapping is reduced by ewith respect to the one corresponding to the initial condition.
3.2. Implementable restart schemes with linear convergence 43 (a) Trajectory of zk.(b) Convergence to optimal objective function value. Figure 3.6: In red, Algorithm 8 applied to Example 3.1. The non-restarted FISTA is depicted in blue. The red circles in (a) depict each of the iterates zk. The black crosses depict each of the restart points rj. However, there is a minor difference with the implementation of this restart scheme that does not fit within the general restart scheme shown in Algorithm 6, which is that FISTA is not restarted using the last point zkthat it generated, but instead the last ykis used. Due to this small difference, we now present the complete pseudocode of this restart scheme. For an initial starting point r0∈Rnz, and an exit tolerance ∈R>0, Algorithm 8 implements FISTA with the restart scheme (3.20) for solving the class of optimization problems (3.9) under Assumtions 2.1 and 3.9. Steps 3 to 6 are the steps of the non-restarted FISTA algorithm (see Algorithm 4), whereas steps 8 to 11 implement the restart procedure and step 7 checks the restart condition (3.20). Notice that in step 9 the new restart point rjis not assigned zk, as would be the case of the general restart scheme described by Algorithm 6. Instead, the current iterate ykis used, both as the new initial condition (see step 9) as for the new “initial” value of the gradient mapping (see step 10). The algorithm generates a sequence ρj=kG(rj)kR−1and exits when ρj≤, i.e., when an -accurate solution of problem (3.9) has been found. This solution is returned in zout, alongside the total number of restart iterations jout and total number of FISTA iterations kout. Example 3.17. Figure 3.6 shows the use of Algorithm 8 to solve Example 3.1 along with the result of the non-restarted FISTA algorithm. The total number of FISTA iterations is kout = 431 and the number of restarts is jout = 14. The convergence properties of this restart scheme are derived from the following proposition, which shows that, under the assumption of the quadratic functional growth condition, kG(y)kR−1can be used to bound the distance of T(y) to the optimal set Ωfof problem (3.9). A similar result can be found in [46, Theorem 7]. However, our proposition can be applied to any non-necessarily feasible point y∈Rnz, i.e., yneed not belong to Z. Analogous results can also be found in other publications, such as [56, 57].
44 Chapter 3. Restart schemes for accelerated first order methods Proposition 3.18. Consider problem (3.9) under Assumtions 2.1 and 3.9 and let Ωfbe the optimal set of problem (3.9). Then, kTf,Z R(y)−πR Ωf(Tf,Z R(y))kR≤1 + √1 + µ µkGf,Z R(y)kR−1,∀y∈Rnz. Proof: To simplify the notation, let us define y+. =Tf,Z R(y), ¯y. =πR Ωf(y) and ¯y+. =πR Ωf(Tf,Z R(y)). From equation (2.7a) in Proposition 2.5 we have that f(y+)−f(z)≤ hG(y), y+−zi+1 2kG(y)k2 R−1,∀y∈Rnz,∀z∈ Z. Particularizing this inequality to z= ¯y+we obtain f(y+)−f(¯y+)≤ hG(y), y+−¯y+i+1 2kG(y)k2 R−1. Since ¯y+∈Ωfwe have f(¯y+) = f∗. Therefore, we infer from Assumtion 3.9 that µ 2ky+−¯y+k2 R≤f(y+)−f∗=f(y+)−f(¯y+) ≤ hG(y), y+−¯y+i+1 2kG(y)k2 R−1 =hR−1 2G(y), R1 2(y+−¯y+)i+1 2kG(y)k2 R−1 (∗) ≤ kG(y)kR−1ky+−¯y+kR+1 2kG(y)k2 R−1, where in (∗) we are making use of the Cauchy-Schwarz inequality. Then, adding 1 2ky+−¯y+k2 Rto both terms of the last inequality, we obtain 1 2(µ+ 1) ky+−¯y+k2 R≤1 2kG(y)kR−1+k¯y+−¯ykR2. That is, pµ+ 1ky+−¯y+kR≤ kG(y)kR−1+k¯y+−¯ykR, (pµ+ 1 −1)ky+−¯y+kR≤ kG(y)kR−1, from where we conclude that ky+−¯y+kR≤1 √1 + µ−1kG(y)kR−1=1 + √µ+ 1 µkG(y)kR−1.
3.2. Implementable restart schemes with linear convergence 45 The following theorem gathers the main convergence properties of the restart scheme presented in Algorithm 8. Claim (i) states the sublinear convergence of {kG(yk)kR−1}≥0under the quadratic functional growth assumption without taking into account the restart procedure; claim (ii) provides an upper bound to the number of iterations required to satisfy the restart condition (3.20) (step 7 of Algorithm 8); and claim (iii) states the linear convergence of the composite gradient mapping under the restart procedure. Finally, claim (iii) provides the iteration complexity. Remark 3.19. We note that we use kG(rj)kR−1≤as the exit condition of Algorithm 8 to simplify the proofs and expressions of Theorem 3.20. In a practical setting, however, the objective is to find an -accurate solution in as fewer iterations as possible. Therefore, the exit condition kG(yk−1)kR−1≤may be used instead (see Remark 2.10). Theorem 3.20 (Convergence properties of Algorithm 8).Consider problem (3.9) under Assumtions 2.1 and 3.9. Then, the sequences generated by Algorithm 8 satisfy: (i) kGf,Z R(yk)kR−1≤4(1 + √µ+ 1) µk kGf,Z R(r0)kR−1,∀k≥1. (ii) kGf,Z R(yk)kR−1≤1 ekGf,Z R(r0)kR−1,∀k≥4e µ(1 + √µ+ 1). (iii) Let a(µ). = max 1,4e(1 + √µ+ 1) µ. Then, min i=0,...,k {kG(yi)kR−1,kG(r0)kR−1} ≤ e1 ejk a(µ)kkG(r0)kR−1. (iv) The total number of iterations kout required to exit the algorithm is upper bounded by kout ≤max 1,4e(1 + √µ+ 1) µlog kG(r0)kR−1 . Proof: Notice that, from step 1 of Algorithm 8, we have that z0=y0=T(r0), which along with Proposition 3.18 implies that kz0−πR Ωf(z0)kR≤1 + √µ+ 1 µkG(r0)kR−1. From Proposition 2.12.(ii) we also have that kG(yk)kR−1≤4 k+ 2kz0−πR Ωf(z0)kR,∀k≥0.
52 Chapter 3. Restart schemes for accelerated first order methods Algorithm 10: Restart scheme for AFOMs based on Ad Require: r0∈dom(f), ∈R>0 1m0←1, m−1←1, j← −1, k←0 2repeat 3j←j+ 1 4sj← sf(rj−1)−f(rj) f(rj−2)−f(rj)if j≥2 0 otherwise 5nj←max{mj,4sjmj−1} 6[rj+1, mj+1]← Ad(rj, nj) 7k←k+mj+1 8until f(rj)−f(rj+1)≤ Output: rout ←rj+1,jout ←j,kout ←k Algorithm 10 shows the proposed restart scheme. Starting from an initial condition r0∈dom(f) and given a AFOM Asatisfying Assumtion 3.24, Algorithm 10 makes successive calls to Ad(Algorithm 9), with a minimum number of iterations njthat is adapted at each restart iteration jto take into account the evolution of the previous objective function values. The algorithm generates a sequence {rj}≥0that converges to an optimal solution z∗of problem (3.22) under Assumtion 3.21 as j→+∞. The algorithm returns an -accurate solution rout in terms of the exit condition shown in step 8 (see Remark 3.32 for some discussion on this exit condition); the total number of restart iterations jout; and the total number of iterations kout performed by algorithm A. The idea behind the restart scheme proposed in Algorithm 10 bares some similarities to the one presented in Section 3.2.1 (see Algorithm 7), in that it makes use of a minimum number of iterations nj, whose value at each iteration j is determined by the evolution of the past objective function values. The similarity is fairly obvious if step 4 of Algorithm 10 is compared with step 7 of Algorithm 7. In both algorithms, the minimum number of iterations is increased if the fraction f(rj−1)−f(rj) f(rj−2)−f(rj−1)(3.26) is larger that a certain amount. In this case, however, the reason behind using (3.26) as a measure for determining the (possible) increase of njdoes not have such an intuitive explanation as the one discussed in Section 3.2.1. Example 3.30. Figure 3.8 shows the use of Algorithm 10 to solve Example 3.1 along with the result of the non-restarted FISTA algorithm. We take Aas MFISTA, which can be implemented as described in Remark 3.28. The total number of MFISTA iterations is kout = 239 and the number of restarts is jout = 5.
3.2. Implementable restart schemes with linear convergence 53 (a) Trajectory of zk.(b) Convergence to optimal objective function value. Figure 3.8: In red, Algorithm 10 applied to Example 3.1 using MFISTA. The non-restarted FISTA is depicted in blue. The red circles in (a) depict each of the iterates zk. The black crosses depict each of the restart points rj. The following proposition provides some key results on the iterates of Algorithm 10. Proposition 3.31. Let Assumtion 3.24 hold and consider Algorithm 10 for a given initial condition r0∈Vf(ρ)⊆dom(f) and an exit tolerance ∈R>0. Then: (i) Proposition 3.29 can be applied to the iterates of Algorithm 10, i.e., taking z0≡rj,n≡nj,zm≡rj+1 and m≡mj+1. (ii) The sequence {mj}≥0is non-decreasing. In particular, mj≤nj≤mj+1,∀j∈Zjout 0.(3.27) (iii) The sequence {sj}≥0satisfies sj∈(0,1], ∀j∈Zjout 2. Proof: Since r0∈Vf(ρ) for some ρ∈R>0, and each rj, for j≥1, is obtained from a call to Algorithm 9 (see step 6), we have in view of Proposition 3.29.(i) that the iterates zkalso satisfy zj∈Vf(ρ), ∀j∈Zjout 0, for the same value of ρ. Therefore, Proposition 3.29 can be applied to each call to Ad, thus proving claim (i). That is, for every j∈Z≥0, the iterates of Algorithm 10 satisfy f(rj+1)≤f(rj)−1 2Lfkg(rj)k2 ∗,(3.28a) f(rj+1)−f∗≤nρ mj+1 + 12 (f(rj)−f∗),(3.28b) nj∈(0,d4nρe]⇒mj+1 ∈[nj,d4nρe].(3.28c) Next, due to step 5 we have mj≤nj,j∈Zjout 0. Moreover, from (3.28c), we have that nj≤mj+1,∀j∈Zjout 0, which proves claim (ii).
54 Chapter 3. Restart schemes for accelerated first order methods Finally, from the exit condition of the algorithm (step 8), we have f(rj−1)−f(rj)> , ∀j∈Zjout 1.(3.29) Additionally, from (3.28a) we have f(rj−2)≥f(rj−1), ∀j∈Zjout 2. Thus, f(rj−2)−f(rj)≥f(rj−1)−f(rj)(3.29) > > 0,∀j∈Zjout 2. Therefore, from step 4, and taking j≥2, we have 0< sj=sf(rj−1)−f(rj) f(rj−2)−f(rj)≤1,∀j∈Zjout 2, which proves claim (iii). Remark 3.32. From Proposition 3.31.(i), we have that we can rearrange Proposition 3.29.(i) to read as kg(rj)k2 ∗≤2Lf(f(rj)−f(rj+1)). Therefore, the exit condition f(rj)−f(rj+1)≤implies kg(rj)k2 ∗≤2Lf. Since, as per Assumtion 3.24, g(rj)serves to characterize the optimality of rj, we conclude that the exit condition of Algorithm 10 also serves to characterize the optimality of rj+1. This means that the exit condition of Algorithm 10 could be replaced by kg(rj)k∗≤˜, where ˜∈Rp. The following theorem presents the main convergence properties of Algorithm 10. In particular, claim (iii) provides its iteration complexity, which we note is of the same order than the optimal iteration complexity (3.2) that can be obtained using AFOMs applied to the class of optimization problems (3.22) under Assumtion 3.21. Theorem 3.33. Let Assumtion 3.24 hold and consider Algorithm 10 for a given initial condition r0∈Vf(ρ)⊆dom(f) and an exit tolerance ∈R>0. Then: (i) The number of calls to Ad(step 6) is bounded. That is, jout is finite. (ii) The number of iterations of Aat each call of Ad(step 6) is upper bounded by d4nρe. That is, mj+1 ≤ d4nρe,∀j∈Zjout 0.(3.30) (iii) The total number of iterations of Arequired to attain the exit condition f(rjout )−f(rjout+1)≤(see step 8), is upper bounded by kout = jout X j=0 mj+1 ≤ed4nρe 25 + 1 log 15 log 1 + f(r0)−f∗ .
3.2. Implementable restart schemes with linear convergence 55 Proof: In the following, we will make use of two technical lemmas: Lemmas 3.35 and 3.36, which we have included, for convenience, immediately after this proof. Let T∈Zbe such that f(rj)−f(rj+1)> , ∀j∈ZT 0,(3.31) is satisfied. Then, defining dj. =f(rj)−f(rj+1), we have f(r0)−f(rT+1) = T X j=0 dj≥(T+ 1) min j=0,...,T dj>(T+ 1). Thus, T+ 1 <f(r0)−f(rT+1) ≤f(r0)−f∗ ≤ρ , from where we infer that the largest integer Tsatisfying (3.31) is bounded. Consequently, the exit condition of Algorithm 10 (step 8) is satisfied within a finite number of iterations, thus proving claim (i). To prove claim (ii), we start by noting that both m1and m2are no larger than d4nρe. Indeed, from step 4 we have that s0=s1= 0, which, in virtue of step 5, implies that n0=m0= 1 and n1=m1. Since n0= 1 is no larger than d4nρe we have from (3.28c) that m1is also upper-bounded by d4nρe. Moreover, since n1=m1≤ d4nρe, we obtain by the same reasoning that m2≤ d4nρe. We now prove that if j≥2 and mj≤ d4nρe, then mj+1 ≤ d4nρe. From step 4 we have s2 j=f(rj−1)−f(rj) f(rj−2)−f(rj)= 1 −f(rj−2)−f(rj−1) f(rj−2)−f(rj) ≤1−f(rj−2)−f(rj−1) f(rj−2)−f∗=f(rj−1)−f∗ f(rj−2)−f∗ (3.28b) ≤nρ mj−1+ 12 . Thus, we have sjmj−1≤nρ. Therefore, nj= max{mj,4sjmj−1} ≤ max{d4nρe,4nρ}=d4nρe, which, along with (3.28c), leads to mj+1 ≤ d4nρe, thus proving the claim. Finally, to prove claim (iii), we start by noting that the computation of each rj+1 is obtained from mj+1 iterations of A(see step 6). Thus, kout = jout X j=0 mj+1 (3.30) ≤(1 + jout)d4nρe.(3.32) Let us denote D. =5 + 1 log 15 log 1 + f(r0)−f∗ . Consider first the case jout < D. Since both jout and Dare integers we infer from this inequality that 1 + jout ≤D. This, along with (3.32), implies that kout ≤ d4nρeD.
56 Chapter 3. Restart schemes for accelerated first order methods Suppose now that jout ≥D. We first recall that Property 3.31.(ii) states that the sequence {mj+1}≥0is non-decreasing. We now rewrite jout as jout =d+tD, where d∈Z0,D−1and t∈Z≥0. Thus, kout = d X j=0 mj+1 + tD X j=1 md+j+1 ≤Dmd+1 +D t X i=1 md+1+iD =D t X i=0 md+1+iD. From Lemma 3.36.(v) we have md+1+iD ≤md+1+(i+1)D √15 ,∀i∈Zt−1 0. Thus, kout ≤D t X i=0 md+1+tD 1 √15t−i . Using now m1+d+tD ≤ d4nρe, see (3.30), we obtain NA Dd4nρe≤ t X i=0 1 √15t−i = t X j=0 1 √15j ≤∞ X j=0 1 √15j =√15 √15 −1≤e 2. Thus, NA≤e 2d4nρeD. Remark 3.34. The iteration complexity provided in Theorem 3.33.(iii) considers the exit condition shown in step 8 of Algorithm 10, whereas the iteration complexities shown in previous sections consider exit conditions of the form kg(rout)k∗≤˜. However, in view of Remark 3.32, this exit condition can also be used in Algorithm 10, in which case the iteration complexity shown in Theorem 3.33.(iii) would be the same but replacing with ˜/(2Lf). The following lemma is used exclusively in the proof of Lemma 3.36. Lemma 3.35. The function ϕ:R→Rdefined as ϕ(s). =1 s2−1max 1,(4s)4, satisfies ϕ(s)≥15,∀s∈(0,√15 4]. Proof: We have that ϕ(s) = 44(s2−s4) if s > 1 4, 1 s2−1 if s≤1 4.
3.2. Implementable restart schemes with linear convergence 57 It is clear that ϕ(·) is monotonically decreasing in (0,1 4]. Thus, min s∈(0,√15 4] ϕ(s) = min s∈[1 4,√15 4] ϕ(s) = min s∈[1 4,√15 4] 44(s2−s4). We notice that the derivative of s2−s4is 2s(1−2s2), which vanishes only once in the interval of interest (at s=1 √2). From here we infer that s2−s4is increasing in [1 4,1 √2) and decreasing in ( 1 √2,√15 4]. Thus, the minimum is attained at the extremes of the interval [1 4,√15 4]. That is, we conclude that min s∈(0,√15 4] ϕ(s) = min{ϕ(1 4), ϕ(√15 4)}= min{15,15}= 15. Lemma 3.36 (A few technical results on the iterates of Algorithm 10).Let Assumtion 3.24 hold and consider Algorithm 10 for a given initial condition r0∈Vf(ρ)⊆dom(f) and an exit tolerance ∈R>0. Assume that jout ≥2 and that there exist T∈Zjout 2and `∈Zjout−T 0such that m`+1 >1 √15m`+1+T. Then, (i) sj∈0,√15 4i,∀j∈Z`+T `+2 , (ii) `+T P j=`+2 log max {1,(4sj)4}<4 log 15, (iii) `+T P j=`+2 log 1 s2 j−1!≤log 1 + f(r0)−f∗ , (iv) T < 5 + 1 log 15 log 1 + f(r0)−f∗ . Additionally, let D. =5 + 1 log 15 log 1 + f(r0)−f∗ . Then, (v) m`+1 ≤1 √15m`+1+D,∀`∈Zjout−D 0.
58 Chapter 3. Restart schemes for accelerated first order methods Proof: From step 4 of Algorithm 10 we have s2 j=f(rj−1)−f(rj) f(rj−2)−f(rj), j ∈Zjout 2. The inequality sj>0, ∀j∈Z`+T `+2 follows from Property 3.31.(iii). In order to prove the first claim it remains to prove the inequality sj≤√15 4,∀j∈Z`+T `+2 . We proceed by reductio ad absurdum. Assume that there is j∈Z`+T `+2 such that sj>√15 4. In this case, mj+1 (3.27) ≥nj= max{mj,4sjmj−1} ≥ 4sjmj−1>√15mj−1. From this and the non-decreasing nature of the sequence {mj}≥0stated in Property 3.31.(ii) we obtain m`+1+T≥mj+1 >√15mj−1≥√15m`+1, which contradicts the initial assumptions, thus proving claim (i). From the non-decreasing nature of the sequence {mj}≥0stated in Property 3.31.(ii) we have mj+1 (3.27) ≥nj= max {mj,4sjmj−1} ≥ mj−1max {1,4sj},∀j∈Z`+T `+2 . Equivalently, log (max {1,4sj})≤ln mj+1 mj−1 ,∀j∈Z`+T `+2 , which implies `+T X j=`+2 log (max {1,4sj})≤ `+T X j=`+2 log mj+1 mj−1 = log m`+Tm`+1+T m`+1m`+2 ≤log m2 `+1+T m2 `+1 = 2 log m`+1+T m`+1 <2 log √15 = log 15. Claim (ii) immediately follows from multiplying this inequality by 4. To prove the third claim, we start by noticing that `+T Y j=`+2 (1 s2 j−1) = `+T Y j=`+2 f(rj−2)−f(rj−1) f(rj−1)−f(rj)=f(r`)−f(r`+1) f(r`+T−1)−f(r`+T). Since `+T≤jout we have f(r`+T−1)−f(r`+T)>>0, which leads to `+T Y j=`+2 (1 s2 j−1) <f(r`)−f(`+1) (3.28a) ≤f(r0)−f(r`+1) ≤f(r0)−f∗ ,
3.3. Numerical results 59 from where claim (iii) directly follows. Next, we add the inequalities given in claims (ii) and (iii) to obtain `+T X j=`+2 log 1 s2 j−1!max 1,(4sj)4!<log 1 + f(r0)−f∗ +4 log 15.(3.33) From claim (i) we have sj∈0,√15 4i,∀j∈Z`+T `+2 . Thus, making use of Lemma 3.35, the left term of (3.33) can be lower bounded by means of the following inequality 15 ≤1 s2−1max 1,(4s)4,∀s∈0,√15 4i. That is, `+T X j=`+2 log 15 <log 1 + f(r0)−f∗ + 4 log 15. Equivalently, (T−1) log 15 <log 1 + f(r0)−f∗ + 4 log 15, thus proving claim (iv). Finally, we prove claim (v) by reductio ad absurdum. If there exist `∈Zjout−D 0 such that m`+1 >1 √15m`+1+D, then we obtain from claim (iv) that D < 5 + 1 ln 15 ln 1 + f(z0)−f∗ , which contradicts the definition of D. 3.3 Numerical results This section presents numerical results comparing the three restart schemes proposed in Section 3.2 with some of the restart schemes of the literature described in Section 3.1. In particular, we compare the restart schemes presented in Algorithms 7, 8 and 10 with the objective function value restart scheme [45], whose restart condition Efis given by (3.6); the gradient alignment restart scheme [45], whose restart condition Egis given by (3.7); and the optimal fixed-rate restart scheme from [46, § 5.2.2], using the restart condition E∗ fgiven by (3.5). For convenience and space considerations, the tables and figures of this section will use the nomenclature of the exit conditions (Ef,Egand E∗ f) to refer to the restart schemes of the literature. Similarly, we will refer to the restart schemes we propose in Section 3.2 by their algorithms (Alg. 7, Alg. 8 and Alg. 10). The results shown here use FISTA (Algorithm 4), with the exception of Algorithm 10, which uses MFISTA (Algorithm 5). We chose MFISTA because, due
60 Chapter 3. Restart schemes for accelerated first order methods to its monotone behavior, step 4 of Algorithm 9 always evaluates to the upper expression (see Remark 3.28). Therefore, MFISTA can be directly used in place of Adin step 6 of Algorithm 10 making use of the restart condition shown in Remark 3.28. The norm of the restart scheme presented in Algorithm 10 will use the same norm k · kRas the other restart schemes. In order to provide a fair comparison between the different schemes, we exit them as soon as an iterate zksatisfying kG(zk)kR−1≤is attained for the selected value of ∈R>0. This way, we can compare how quickly each restart scheme finds an -accurate solution of the problem at hand. Restart scheme E∗ frequires knowledge of f∗, which we obtain by solving the optimization problem using Algorithm 7 with and exit tolerance of = 10−8. As additional numerical results, Section 5.8.7 shows the application of the above mentioned restart schemes to solve the optimization problems of different MPC formulations using FISTA. 3.3.1 Application to Lasso problems This section presents the result of applying the restart schemes to weighted Lasso problems of the form min z∈Rnz 1 2NkAz −bk2 2+kWzk1,(3.34) where z∈Rnz,A∈RN×nzis sparse with an average of 90% of its entries being zero (sparsity is generated by setting a 0.9 probability for each element of the matrix to be 0), nz> N, and b∈RN. Each nonzero element in Aand bis obtained from a Gaussian distribution with zero mean and variance 1. W∈Dnz ++ is a diagonal matrix with elements obtained from a uniform distribution on the interval (0, α]. We note that problems (3.34) can be reformulated in such a way that they satisfy the quadratic growth condition [46, § 6.3]. We take Ras the diagonal matrix constructed as R(i,i)= n X j=1 |H(i,j)|,(3.35) where H=1 NA>A, which due to the Gershgorin Circle Theorem [60, § 7.2] (see also [27, § 6]) satisfies the smoothness condition given in Assumtion 2.1.(ii). Table 3.1 shows the results of solving 100 randomly generated problems (3.34) that share the values of N= 600, nz= 800, and α= 0.003. We take = 10−7. Figure 3.9 and Figure 3.10 show the evolution of kG(zk)kR−1and f(zk)−f∗, respectively, of each one of the restart schemes for one of the Lasso problems (3.34) used to obtain the results of Table 3.1. Additionally, they show the result of applying FISTA without a restart scheme. Table 3.2 and Figures 3.11 and 3.12 show analogous results to Table 3.1 and Figures 3.9 and 3.10, respectively, but taking N= 100, nz= 200 and α= 0.3.
3.3. Numerical results 61 Restart scheme Alg. 7 Alg. 8 Alg. 10 EfEgE∗ f Avg. Iter. 914.69 1431.2 897.46 946.62 894.21 1350.3 Med. Iter. 907 1410.5 880.5 945.5 880.5 1331 Max. Iter. 1155 1800 1131 1349 1199 1696 Min. Iter. 704 1154 737 692 702 1084 Table 3.1: Comparison between restart schemes applied to FISTA to solve 100 problems (3.34) with N= 600, nz= 800, α= 0.003 and = 10−7. 0 320 640 960 1280 1600 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 Figure 3.9: Evolution of the dual norm of the composite gradient mapping for different restart schemes applied to a randomly generated problem (3.34) with N= 600, nz= 800, α= 0.003 and = 10−7. Black dots represent the iterations in which a restart occurred. 0 320 640 960 1280 1600 10 -11 10 -9 10 -7 10 -5 10 -3 10 -1 Figure 3.10: Evolution of the distance, in terms of function value, to the optimal solution for different restart schemes applied to a randomly generated problem (3.34) with N= 600, nz= 800, α= 0.003 and = 10−7. Black dots represent the iterations in which a restart occurred.
69 Chapter 4 Preliminaries This chapter addresses some of the preliminary concepts and ideas related to model predictive control (MPC) and its implementation in embedded systems. We start by introducing the problem formulation in Section 4.1, i.e., by presenting the class of system and system model under consideration as well as the control objective. We then provide a brief description of MPC in Section 4.2 so as to introduce the concepts and nomenclature used throughout the remainder of the dissertation. Finally, in Section 4.3 we present a brief (and non-exhaustive) review of the literature regarding the implementation of MPC in embedded systems. Before we proceed, we provide for the sake of clarity our definition of an embedded system. Some examples of well known embedded systems are Programmable Logic Controllers (PLC), FPGAs, Arduino and Raspberry Pi. Definition 4.1 (Embedded system).We understand by embedded system any computer system that combines both hardware and software and that has a dedicated function within a larger system. In particular, we circumscribe our definition to those electronic systems that are based either on microcontrollers or microprocessors and that are used to control the physical device they are embedded in. Integral to our definition is the fact that embedded systems are characterized by their low computational and memory resources when compared to more powerful computing devices, such as the average desktop PC. 4.1 Problem formulation We consider a discrete-time, linear, time-invariant, state space model x(t+ 1) = Ax(t) + Bu(t) (4.1) where x(t)∈Rnand u(t)∈Rmare the state and control input of the system at the discrete time instant t∈Z,A∈Rn×n, and B∈Rn×m. Additionally, we assume that the system is subject to state and input constraints (x(t), u(t)) ∈ Y ⊆ Rn×Rm,∀t, (4.2)
70 Chapter 4. Preliminaries where the set Ywill be given either by box constraints of the form Y={(x, u)∈Rn×Rm:x≤x≤x, u ≤u≤u},(4.3) where x∈Rn,x∈Rn,u∈Rm, and u∈Rm; or by coupled input-state constraints Y={(x, u)∈Rn×Rm:y≤Ex +Fu ≤y},(4.4) where y∈Rp,y∈Rp,E∈Rp×nand F∈Rp×m. We make the following assumption on model (4.1) and its constraints. Assumption 4.2. We assume that: (i) Model (4.1) is controllable. (ii) The constraint set Yhas a non-empty interior. Remark 4.3. If the constraints are given by (4.3), then Assumtion 4.2.(ii) implies that x< x and u < u. The control objective is to steer the system (4.1) to the given reference (xr, ur)∈ Rn×Rmwhilst satisfying the system constraints (4.2), where we assume that the reference is a (piecewise) constant set-point, i.e. a constant reference whose value can change at any unpredictable time. Obviously, this objective is only possible if the reference is an admissible steady state (see the following definition) of system (4.1) under (4.2). Definition 4.4 (Admissible steady state).Consider a system (4.1) subject to (4.2) under Assumtion 4.2. We say that a pair (ˆx, ˆu)∈Rn×Rmis an admissible steady state of the (4.1) under (4.2) if ˆx=Aˆx+Bˆu(i.e., it is a steady state) and (ˆx, ˆu)∈ Y. We say that it is strictly admissible if (ˆx, ˆu)∈int(Y). If the reference is not an admissible steady state, then we wish to steer the system to the “closest” admissible steady state, for some criterion of closeness. In this dissertation we are concerned with linear systems (4.1). However, for future purposes, we now introduce some notation for the real system from where the linear model is derived. In general, model (4.1) will be a linear representation of some real system governed by a continuous-time, non-linear, ordinary differential equation dx dt=f(x(t),u(t)),(4.5) where x(t)∈Rnand u(t)∈Rmare the state and control input of the system at the continuous time instant t∈R, respectively, and f:Rn×Rm→Rnis a continuous differentiable function. Additionally, the constraints (4.2) will have been obtained so as to satisfy some real system constraints (x(t),u(t)) ∈Y⊆Rn×Rm,∀t. (4.6)
4.2. A brief introduction to model predictive control 71 In the case of box constraints, the set Ywill be given by Y={(x,u)∈Rn×Rm:x≤x≤x,u≤u≤u}. where x∈Rn,x∈Rn,u∈Rm, and u∈Rm. In particular, we consider that model (4.1) has been obtained by linearizing the real system model (4.5) around an operating point (x◦,u◦), which we assume to be an admissible steady state of (4.5) under (4.6). Additionally we consider that xand umay have been subject to scaling using the scaling matrices Nx∈Dn and Nu∈Dmas follows: x=Nx(x−x◦), u =Nu(u−u◦). There are many other ways in which a linear model can be derived from a non-linear one or from data of the system [62]. Our above consideration is merely the one we take in the numerical examples we show in this dissertation. We refer to x,uand their bounds as being in incremental units to emphasize that they are (scaled) increments with respect to the operating point of the real system model. Additionally, we refer to x,uand their bounds as being in engineering units, and we will represent them using calligraphic font. We use this naming convention to emphasize the fact that they relate to the real system (even though, in many occasions, the state has no tangible physical interpretation). 4.2 A brief introduction to model predictive control Model predictive control (MPC) is an advanced control strategy that is very prevalent in the current control literature due to its inherent ability to provide constraint satisfaction and asymptotic stability to the given target [63, 64]. In MPC, the control action is obtained, at each sample time, from the solution of an optimization problem in which a prediction model is used to predict the future evolution of the system over a given prediction horizon. The need for solving an optimization problem in real-time and the fact that a model of the system is required have historically hindered the use of MPC in many industries and applications. In spite of this, there are numerous publications on the successful use of MPC in different areas, including water distribution networks [65], micro grids [66], power converters [67], automotive applications [68, 69], or temperature control of buildings [70], data centers [71] and of the heat, ventilation and air conditioning systems [17, 72, 73]. Different MPC formulations and control architectures are used in the above references. In this dissertation we focus on linear MPC formulations with simple constraints. This section is not intended as a thorough description nor explanation of MPC. Instead, it is a brief introduction to the class of MPC problems that we focus on in this dissertation, as well as to the nomenclature that we will use.
72 Chapter 4. Preliminaries We consider linear MPC formulations described by the following parametric optimization problem: min x,u J. = N−1 X j=0 `j(x,u;xr, ur) + Vt(x,u;xr, ur) (4.7a) s.t. xj+1 =Axj+Buj, j ∈ZN−1 0(4.7b) x0=x(t) (4.7c) (xj, uj)∈ Yj, j ∈ZN−1 0(4.7d) xN∈ Xt,(4.7e) where N∈Z>0is the prediction horizon; (4.7b) is the prediction model (4.1), with A∈Rn×nand B∈Rn×m; the sequences x= (x0, . . . , xN) and u= (u0, . . . , uN−1) are the predicted states xj∈Rnand inputs uj∈Rmthroughout the prediction horizon; xr∈Rnand ur∈Rmare the given state and input reference; x(t)∈Rn is the system state at the current discrete time instant t;Yj⊆Rn×Rmare the constraints for the prediction step j;`j(·) and Vt(·) are the stage costs for each prediction step jand the terminal cost functions, respectively, which in general are functions of the predicted states xand inputs uand are parametrized by the reference (xr, ur); and Xtis the terminal set. Constraint (4.7e) is known as the terminal constraint. The MPC formulations that we discuss and present in the following chapters loosely fall within the above general formulation, although some of them include additional decision variables and constraints. However, we always consider that the functions `jand Vtof the cost function J(4.7a) are real valued convex functions and that Xt⊆Rnis a closed convex set. Therefore, problem (4.7) is a convex optimization problem. We denote an optimal solution of problem (4.7) (assuming that one exists) by x∗= (x∗ 0, x∗ 1, . . . , x∗ N) and u∗= (u∗ 0, u∗ 1, . . . , u∗ N−1). The MPC control law, i.e., the control action u(t) to be applied to the system each discrete time instant t, is given by u(t) = u∗ 0. At the next sample time t+1, the MPC optimization problem is recalculated taking the new system state x(t+ 1) in (4.7c), and so forth. The following two definitions play an important role in the stability of MPC. Definition 4.5 (Feasibility region).Consider an MPC formulation defined by some optimization problem parametrized by the current system state at time instant t∈Z. Its feasibility region is the set of states x(t)∈Rnfor which said optimization problem is feasible. Definition 4.6 (Domain of attraction).Consider an MPC formulation parametrized by the current system state at time instant t∈Zand whose objective is to steer the system to a reference xr∈Rn. Its domain of attraction is the set of states x(t)∈Rnsuch that the MPC controller steers the closed-loop system to the reference xrwhile satisfying the constraints.
4.3. A brief state of the art 73 We will often consider the following standard assumption. Assumption 4.7. Let Assumtion 4.2 hold and assume that (i) (xr, ur) is an admissible steady state of system (4.1) under (4.2). (ii) x(t) strictly belongs to the feasibility region (Definition 4.5) of the MPC formulation, i.e., there exist xand usatisfying (4.7b) and (4.7c) for which (4.7d) and (4.7e) are strictly satisfied. In MPC, the terminal set Xtis typically taken as an admissible invariant set (see Definition N.8) of the system (4.1) for the constraints (4.2) under some terminal control law. That is, for some predefined control law, any state belonging to Xtwill remain in Xt. The use of a terminal admissible invariant set provides the MPC controller with stability guarantees if its other ingredients are suitably designed [64]. The control objective of the MPC controller is to steer the system to the reference (xr, ur) whilst satisfying the system constraints (4.2), which will only be possible if (xr, ur) is an admissible steady state of (4.1) under (4.2). Remark 4.8. We note that our control objective, as stated above, is to steer the system to an admissible steady state of the linear model (4.1) whilst satisfying (4.2). Since this model is a linearization of a (possibly) non-linear model, our control objective does not guarantee that the real system will converge to the desired set point, unless additional ingredients are added to the overall control architecture [8, 74]. During this dissertation we will not delve into this topic any further, since our main interest is in the development of MPC formulations and solvers suitable for their implementation in embedded systems, and not on the additional ingredients required for their use in a real environment. 4.3 A brief state of the art The implementation of linear MPC in embedded systems is an extensively researched topic in the field of control, with many publications tackling different embedded systems, MPC formulations and approaches. This section presents a (non-exhaustive) review of the state of the art. We focus on the implementation of linear MPC, instead of non-linear, because the solvers proposed in this dissertation all fall into this paradigm. There are also multiple results on the implementation of nonlinear MPC in embedded systems [75, 76, 77, 78], as well as multiple tools and solvers [79, 80]. Additionally, we focus on the implementation of nominal MPC, although results of the implementation of robust MPC are also available, including the min-max robust MPC approach [81, 82, 83]. Henceforth, we will refer to “linear MPC” simply as “MPC”, unless specifically stated otherwise. The implementation of MPC in embedded systems can be divided into two main paradigms: explicit MPC, and online MPC.
74 Chapter 4. Preliminaries In explicit MPC [84, 85], the solution of the MPC control law is computed offline and stored in the embedded system. In [86] the authors show that the control law of standard MPC subject to state and input constraints is given by a piece-wise affine and continuous function of the state. This function can be computed offline and stored as a look-up table. Therefore, no optimization problem has to be solved online, since the control action is determined by evaluating the look-up table. This approach has been used in various publications to implement MPC in embedded systems [87, 88, 89, 90, 91]. However, its drawback is that the size of the look-up table, and thus the memory requirements and computational time needed to evaluate it, becomes prohibitively large for medium to large-sized systems and/or for problems involving many constraints. The second main paradigm is to solve the MPC’s optimization problem at each sample time for the current system state. A large portion of the results in this field come from the recent development of solvers that are tailored to their use in embedded systems, which can be used to solve the MPC optimization problem online. Many of them generate library-free code that is tailored to the problem to be solved, thus resulting in a rather efficient implementation. There are a wide range of solvers, each particularized to a certain class of optimization problem and based on a particular optimization method. A few of the most noteworthy and well known solvers that have been used for the embedded implementation of MPC are the following: CVXGEN [92] is an web-based “code generator for convex optimization problems that can be reduced to quadratic programs” [92] and that is better suited for small scale problems. It generates library-free C code tailored to the specific QP problem, which is solved using a primal-dual interior point method. qpOASES [93, 94], is an open-source C++ QP solver based on a parametric active set method that, according to the authors, is “particularly suited for model predictive control applications” [94]. OSQP [35] is open-source sparse operator splitting solver for QP problems that is based on the ADMM algorithm. The solver is written in C, but has interfaces to various other programming languages. ODYS [95] is a proprietary solver for QP problems with a specific focus on MPC. It has been used to implement MPC in various industries, such as, for instance, in automotive mass production [69]. FiOrdOs [96] is a Matlab toolbox for automated C-code generation of first order methods for parametric convex programs that is well suited for its application to MPC. The optimization problem is solved using the gradient method or the fast gradient method.
4.3. A brief state of the art 75 FORCES [97] is a proprietary solver that generates code for convex optimization problems that is based on interior point methods and well suited for MPC. qpDUNES [98] is an C-based open-source QP solver based on a dual Newton strategy. Some examples of these tools being used to implement MPC in industrial embedded systems include [99], where qpOASES is implemented in a PLC to control a MISO system; [100], where qpOASES, CVXGEN and FiOrdOs are implemented in a PLC and compared (see also [101] for another implementation using FiOrdOs); or [102] which compares the scalability of FiOrdOs and qpOASES in a PLC. A very useful tool for building optimizations problems and solving them with the above solvers (and many others) is the YALMIP toolbox for Matlab [103]. The above solvers are not specific to MPC, although some consider and/or provide tools for it. As such, even though they can be used to solve a wider class of optimization problems, their use for MPC has been reported in numerous publications. However, another approach followed by some authors is to develop solvers that are tailored or more focused on MPC. These solvers may attain a more efficient implementation due to their narrow focus. Some noteworthy examples include: µAO-MPC [104] is a code generation tool for MPC subject to box constraints on the inputs. It generates self-contained code tailored to the specific MPC problem and system. The optimization problem is solved using the augmented Lagrangian method along with the fast gradient method. PRESAS [105] is a primal active-set solver for block-sparse QPs with a particular focus on MPC (both linear and non-linear). HPMPC [106] is a C library for the implementation of MPC. Various routines are used, depending on if the MPC is unconstrained or constrained (Riccati recursion, interior point and ADMM). Its successor, HPIPM [107], which is currently under active development, focuses on interior point methods for QPs, with a particular focus on MPC. SPCIES [4] is the Matlab toolbox that contains the solvers that have been developed as part of this dissertation. We now briefly describe other noteworthy publications in this field that don’t quite fall into the above paradigms, either because they implement MPC in embedded systems but without presenting a fully-fledged solver, or because they present novel methods/approaches for its implementation. In [108], an online implementation of embedded MPC is presented, where the MPC optimization problem is solved using a combination of Nesterov’s fast gradient method and method of multipliers. A parallel coordinate descent algorithm
76 Chapter 4. Preliminaries for network systems is presented in [109] and applied to MPC in a PLC. In [110], an accelerated dual gradient projection algorithm for embedded MPC with an iteration complexity that grows linearly with the prediction horizon is presented. The article [111] presents the implementation of MPC in FPGAs using the fast gradient method and ADMM. In [112], Nesterov’s fast gradient method is used to solve MPC in a PLC. A tool for software-hardware code generation of operator splitting methods for MPC in FPGAs is presented in [113], The application of an accelerated dual proximal algorithm is proposed in [114] to finding the solution of the infinite horizon constrained LQR problem. In [115], a primal-dual iterative algorithm particularized to MPC is developed for its implementation in FPGAs. In [116], MPC is programmed on a FPGA using high-level synthesis tools along with µAO-MPC. A primal active-set method for bounded-variable least-squares problems [117] is implemented in a PLC to solve MPC in [15]. An ADMM algorithm for MPC which exploits the problem structure by decoupling the system dynamics is presented in [118]. Decoupling the system dynamics allows for a parallel implementation of the algorithm. In [119], another ADMM based solver for MPC that decouples the system dynamics is presented, in this case by exploiting the symmetry of the system. As evidences from the above references, research in this field is very extensive, encompassing a wide range of embedded systems, optimization methods and approaches. For additional references we refer the reader to the surveys [120, 121, 122, 123].
77 Chapter 5 Sparse solvers for model predictive control This chapter considers system (4.1) subject to the box constraints (4.3). This chapter presents sparse solvers for various MPC formulations in which the choice of formulation and optimization algorithm leads to an efficient solver that exploit the structure of the problem. That is, we define the MPC formulations and the decision variables in such a way that the ingredients of the first order methods (i.e., the resulting matrices) are simple enough to be able to perform matrix-vector operations are simple and can be performed without having to use the typical sparse matrix representations (e.g., the compressed sparse column/row or dictionary of keys formats). The advantage of doing this is that we do not need to store all the non-zero elements of the matrices nor the arrays that determine their positions within them. Instead, we only need to store the repeating submatrices once, and the sparse matrix-vector operations are performed by direct identification of their inner structure. This results in solvers with low iteration complexity and memory footprint; although, in some aspects, at the expense of some restrictions. The solvers presented in this chapter have been included in the SPCIES toolbox for Matlab1[4], which is available at https://github.com/GepocUS/Spcies. We do not provide the exact pseudocode of the solvers, since the specific ways in which we perform the sparse matrix-vector operations, as well as certain details such as what variables are stored, are subject to future changes and improvements. Instead, we describe the pseudocode in more general terms. We show how the sparse matrix-vector operations arise, the inner structure of the matrices, and providing enough information to understand which variables are required and which are not. The reader can find the exact code of the solvers in [4]. 1For future reference, the version of the toolbox as of the date of presentation of this dissertation is v0.3.2. Future versions may have made improvements upon the results shown here.
84 Chapter 5. Sparse solvers for model predictive control For a given λ, let zλ∈Rnzbe defined as zλ= arg min z∈Z L(z, λ) (5.8) It is obvious that φ(λ) = L(zλ, λ). Since Jis a strongly convex function, we have that, under Assumtion 5.11.(iv), the dual function ψis a continuously differentiable smooth concave function [124, Lemma 3.1] whose gradient is given by [19, Example 2.2.4] ∇ψ(λ) = b−Gzλ.(5.9) The following proposition characterizes the smoothness of ψin terms of the matrix W. =GH−1G>(c.f., Assumtion 2.1.(ii)). Its proof is based on the following lemma. Note that, since H∈Sn ++ and rank(G) = mz, we have that W∈Smz ++. Lemma 5.13. Consider the primal and dual problems (5.5) and (5.6), where the dual function ψis given by (5.2.1). Let Assumtion 5.11 hold and zλbe given by (5.8). Assume that ∆z∈Rnzis such that zλ+ ∆z∈ Z. Then, h∆z, Hzλ+q−G>λi ≥ 0. Proof: From the convexity of Z, have that (1−µ)zλ+µ(zλ+∆z)∈ Z,∀µ∈[0,1], (see Definition N.2) which leads to zλ+µ∆z∈ Z,∀µ∈[0,1]. From the optimality of zλ, we have ψ(λ) = 1 2z> λHzλ+q>zλ+hλ, b −Gzλi ≤1 2(zλ+µ∆z)>H(zλ+µ∆z) + hq, zλ+µ∆zi+hλ, b −G(zλ+ ∆z)i =ψ(λ) + µh∆z, Hzλ+q−G>λi+µ2 2∆z>H∆z, ∀µ∈[0,1]. Therefore, µh∆z, Hzλ+q−G>λi+µ2 2∆z>H∆z≥0,∀µ∈[0,1], which is trivially satisfied if µ= 0. Otherwise, dividing by µ, we have that h∆z, Hzλ+q−G>λi+µ 2∆z>H∆z≥0,∀µ∈(0,1]. For this to be true ∀µ∈(0,1], it must be true that h∆z, Hzλ+q−G>λi ≥ 0. Proposition 5.14 (Smoothness of the QP’s dual function).Consider the primal and dual problems (5.5) and (5.6), where the dual function ψis given by (5.2.1). Let Assumtion 5.11 hold and zλbe given by (5.8). Denote W. =GH−1G>. Then for any (λ, ∆λ)∈Rmz×Rmz, ψ(λ+ ∆λ)≥ψ(λ) + h∆λ, b −Gzλi− 1 2k∆λk2 W.
5.2. Solving QPs with FISTA through duality 85 Proof: Let us define ∆zas zλ+ ∆z= arg min z∈Z 1 2z>Hz +q>z+hλ+ ∆λ, b −Gzi. It is clear that zλ+ ∆z∈ Z. Then, from the definition of ψ, we have ψ(λ+∆λ) = 1 2(zλ+∆z)>H(zλ+∆z) + hq, zλ+∆zi+hλ+∆λ, b −G(zλ+∆z)i =ψ(λ) + 1 2∆z>H∆z+h∆z, Hzλ+q−G>λi+h∆λ, b −G(zλ+∆z)i (∗) ≥ψ(λ) + 1 2∆z>H∆z+h∆λ, b −G(zλ+∆z)i =ψ(λ) + h∆λ, b −Gzλi+1 2∆z>H∆z+h∆λ, −G∆zi ≥ψ(λ) + h∆λ, b −Gzλi+ min ∆z1 2∆z>H∆z+h∆λ, −G∆zi (∗∗) =ψ(λ) + h∆λ, b −Gzλi− 1 2∆λ>GH−1G>∆λ, where (∗) is due to Lemma 5.13 and (∗∗) follows from min ∆z1 2∆z>H∆z+h∆λ, −G∆zi=−1 2∆λ>GH−1G>∆λ, which is a well known result in the field of convex optimization. 5.2.2 Solving the QP’s dual problem with FISTA We now show how to solve the dual problem (5.6) using FISTA (Algorithm 4). The results we now present are well known and can be found in several prior publications, but are included here for completeness. Since FISTA is expressed in terms of the minimization of a convex function, instead of maximization of a concave one, we consider the convex problem −min λ∈Rmz−ψ(λ).(5.10) We note that, from (5.9), we have ∇(−ψ(λ)) = −(b−Gzλ).(5.11) Then, under Assumtion 5.11, and taking into consideration Proposition 5.14, we have that Assumtion 2.1 is satisfied for (5.10). In particular, problem (5.10) has no non-smooth term Ψ, and the matrix Rthat characterizes its smoothness (Assumtion 2.1.(ii)) is W. =GH−1G>. Therefore, the iterates of FISTA applied to problem (5.10) will converge to λ∗.
86 Chapter 5. Sparse solvers for model predictive control Algorithm 14: FISTA for solving QP problem (5.5) Require: λ∈Rmz,∈R>0 1t0←1, k←0 2qk←q−G>λ 3zk←arg min z∈Z 1 2z>Hz +q> kz 4Γk← −(Gzk−b) 5∆λk←solution of W∆λ= Γk 6yk←λ+ ∆λk 7λk←λ+ ∆λk 8repeat 9k←k+ 1 10 qk←q−G>yk−1 11 zk←arg min z∈Z 1 2z>Hz +q> kz 12 Γk← −(Gzk−b) 13 ∆λk←solution of W∆λ= Γk 14 λk←∆λk+yk−1 15 tk←1 21 + q1+4t2 k−1 16 yk←λk+tk−1−1 tk (λk−λk−1) 17 until kΓkk∞≤ Output: ˜z∗←zk,˜ λ∗←yk−1 Let us now take a closer look at step 4 of Algorithm 4, which performs the assignment λk← T −ψ,Rmz W(yk−1). That is, λk←arg min λ∈Rmzh∇(−ψ(yk−1)), λ −yk−1i+1 2kλ−yk−1k2 W = arg min λ∈RmzhGzyk−1−b, λ −yk−1i+1 2kλ−yk−1k2 W, which is a strongly convex unconstrained QP problem. Therefore, λkis the solution of the system of equations W(λk−yk−1) = −(Gzyk−1−b). By defining ∆λk=λk−yk−1, step 4 of Algorithm 4 reduces to first computing zyk−1, then solving the system of equations W∆λk=−(Gzyk−1−b), and then performing the assignment λk←∆λk+yk−1. FISTA algorithm applied to the dual problem (5.10) is shown in Algorithm 14. Note that its exit condition (step 17) does not correspond to the exit condition of Algorithm 4, which would be kyk−1−λkkW≤. We use the exit condition kΓkk∞≤because it is computationally cheaper to evaluate, since Γk∈Rmzis
5.3. Solving QPs with ADMM 87 Algorithm 15: ADMM for solving QP problem (5.12) Require: v0∈Rnz,λ0∈Rnz,ρ∈R>0,p∈R>0,d∈R>0 1k←0 2repeat 3qk←q+λk−ρvk 4zk+1 ←min z 1 2z>Hρz+q> kz, s.t. Gz =b 5ˆqk← −ρzk+1 −λk 6vk+1 ←arg min v∈Z ρ 2v>v+ ˆq> kv 7λk+1 ←λk+ρ(zk+1 −vk+1) 8k←k+ 1 9until kzk−vkk∞≤pand kvk−vk−1k∞≤d Output: ˜z∗←zk, ˜v∗←vk,˜ λ∗←λk already computed in step 12. The fact that Γkis a measure of optimality of the current iterate follows from Γk=−(Gzk−b) = −(Gzyk−1−b)(5.11) =−∇(−ψ(yk−1)). Then, since (5.10) is an unconstrained smooth convex problem, we have that λ∗=yk−1⇐⇒ ∇(−ψ(yk−1)) = 0. That is, λ∗=yk−1⇐⇒ Γk= 0. Additionally, from Proposition 5.12, we have that z∗=zλ∗. Thus, the outputs of the algorithm are the suboptimal dual solution ˜ λ∗=yk−1and the suboptimal primal solution ˜z∗=z˜ λ∗=zyk−1=zk, where the suboptimality is determined by . Remark 5.15. Note that steps 3 and 11 of Algorithm 14 require solving a constrained QP problem, which in general is not trivial. However, in the following sections we will consider assumptions under which this problem will have a simple and explicit solution. Additionally, steps 5 and 13 require solving a system of equations, which can be computationally demanding. However, we will consider systems of equations that satisfy Assumtion 5.1, and may therefore be sparsely solved using Algorithm 11. 5.3 Solving QPs with ADMM This section explains how problem (5.5) under Assumtion 5.11 can be solved using ADMM (Algorithm 2). The results we now present are well known and can be found in several prior publications, but are included here for completeness. We start by rewriting (5.5) into a problem of class (2.12) by taking min z,v 1 2z>Hz +q>z(5.12a) s.t. Gz =b(5.12b) v∈ Z (5.12c) z−v= 0.(5.12d)
88 Chapter 5. Sparse solvers for model predictive control That is, the ingredients of (2.12) are given by: C=Inz,D=−Inz,d=0nz, f(z) = 1 2z>Hz +q>z+IGz=b(z), g(v) = IZ(v), where IGz=bis the indicator function of the set {z∈Rnz:Gz =b}. Let us revisit steps 3 and 4 of Algorithm 2 when applied to (5.12). Step 3 of Algorithm 2 now reads as zk+1 ←arg min z 1 2z>Hz +q>z+ρ 2kz−vk+1 ρλkk2 2 s.t. Gz =b = arg min z 1 2z>Hρz+q> kz(5.13) s.t. Gz =b, where Hρ=H+ρInzand qk=q+λk−ρvk. Step 4 of Algorithm 2 now reads as vk+1 ←arg min v∈Z ρ 2kzk+1 −v+1 ρλkk2 2 = arg min v∈Z ρ 2v>v+ ˆq> kv, (5.14) where ˆqk=−ρzk+1 −λk. Algorithm 15 shows the result of particularizing Algorithm 2 to problem (5.12) following the above discussion. Remark 5.16. As in Remark 5.15, steps 4 and 6 of Algorithm 15 require solving an equality-constrained QP and a constrained QP, respectively. In general, these problems are not necessarily simple to solve. However, in this chapter we will consider assumptions under which both steps can be solved using Algorithms 12 and 13, respectively. 5.4 Simple standard MPC formulations This section considers two simple standard MPC formulations, which are particularizations of the general MPC formulation (4.7). We describe the formulation, and present sparse solvers based on the algorithms shown in Sections 5.2 and 5.3. The solvers were originally presented in [8, 9], although with a bigger emphasis on their implementation in PLCs.
5.4. Simple standard MPC formulations 89 The first one is an MPC formulation with a terminal equality constraint: min x,u N−1 X j=0 kxj−xrk2 Q+kuj−urk2 R(5.15a) s.t. xj+1 =Axj+Buj, j ∈ZN−1 0(5.15b) x0=x(t) (5.15c) x≤xj≤x, j ∈ZN−1 1(5.15d) u≤uj≤u, j ∈ZN−1 0(5.15e) xN=xr,(5.15f) and the second one is an MPC formulation without a terminal constraint: min x,u N−1 X j=0 kxj−xrk2 Q+kuj−urk2 R+kxN−xrk2 T(5.16a) s.t. xj+1 =Axj+Buj, j ∈ZN−1 0(5.16b) x0=x(t) (5.16c) x≤xj≤x, j ∈ZN 1(5.16d) u≤uj≤u, j ∈ZN−1 0.(5.16e) We consider Assumtion 4.7 to hold for both the above MPC formulations. Assumptions on the cost function matrices Q,Rand Twill be stated further ahead when applicable. The advantage of the MPC formulation (5.15) is that it avoids the computation of an admissible invariant set Xt(Definition N.8), which can be very demanding even for average-sized systems [125], by using the simplest one: the singleton {xr}To see that under Assumtion 4.7, the singleton {xr}is, indeed, an admissible invariant set of the system, note that the admissible terminal control law u(t) = ur keeps the system at xr. However, this comes at the expense of a reduction of the domain of attraction of the MPC controller (Definition 4.6) due to the use of the smallest admissible invariant set of the system. This drawback can be avoided by simply eliminating the terminal constraint, leading to the MPC formulation (5.16). However, this comes at the expense of losing the stability guarantees that accompany the use of a suitable terminal constraint. In particular, the controller may not stabilize every feasible initial state; only those contained within a certain region which may be difficult to characterize [126]. Thus, the feasibility region and the domain of attraction of the MPC formulation (5.16) may not be the same. Remark 5.17. We note that the MPC formulation (5.15) does not include a terminal cost because it does not require one for stability purposes and we chose to focus on a simple MPC formulation in order to attain solvers with low iteration complexity. The MPC formulations (5.15) and (5.16) can be posed as QPs (5.5) as follows.
90 Chapter 5. Sparse solvers for model predictive control QP problem of the MPC formulation (5.15): Take the decision variables as: z= (u0, x1, u1, x2, u2, . . . , xN−1, uN−1).(5.17a) Then, the resulting ingredients of the QP problem (5.5) are given by H= diag(R, Q, R, Q, R, . . . , Q, R),(5.17b) q=−(Rur, Qxr, Rur, Qxr, Rur, . . . , Qxr, Rur),(5.17c) G= B−In0··· ··· 0 0 0A B −In··· 0 0 0 0 .........0 0 0 0 ··· A B −In0 0 0 ··· 0 0 A B ,(5.17d) b= (−Ax(t),0n,0n,...,0n, xr),(5.17e) Z={z:z≤z≤z},(5.17f) where z= (u, x, u, x, u, . . . , x, u),(5.17g) z= (u, x, u, x, u, . . . , x, u).(5.17h) The dimensions of the QP problem are nz= (N−1)(n+m) + mand mz=Nn. QP problem of the MPC formulation (5.16): Take the decision variables as: z= (u0, x1, u1, x2, u2, . . . , xN−1, uN−1, xN).(5.18a) Then, the resulting ingredients of the QP problem (5.5) are given by H= diag(R, Q, R, Q, R, . . . , Q, R, T),(5.18b) q=−(Rur, Qxr, Rur, Qxr, Rur, . . . , Qxr, Rur, Txr),(5.18c) G= B−In0··· ··· 0 0A B −In··· 0 0 0 .........0 0 0 0 A B −In ,(5.18d) b= (−Ax(t),0n,0n,...,0n),(5.18e) Z={z:z≤z≤z},(5.18f) where z= (u, x, u, x, u, . . . , x, u, x),(5.18g) z= (u, x, u, x, u, . . . , x, u, x).(5.18h) The dimensions of the QP problem are given by nz=N(n+m) and mz=Nn.
5.4. Simple standard MPC formulations 91 5.4.1 FISTA-based solver for standard MPC We apply the results of Section 5.2 to the QP problems of the MPC formulations (5.15) and (5.16), whose ingredients are given by (5.17) and (5.18), respectively, under the following assumption. Assumption 5.18. Let Assumtion 4.7 hold and assume that Q∈Dn ++,R∈Dm ++ and T∈Dn ++. It is easy to see that the satisfaction of Assumtion 5.18 implies the satisfaction of Assumtion 5.11 for both QP problems, where, additionally, H∈Dnz ++. Therefore, Algorithm 14 can be applied to them. Let us take a closer look at the computationally expensive steps of Algorithm 14. Steps 5 and 13 of Algorithm 14 require solving W∆λ= Γk. However, due to the banded structures of Gand H, matrix W=GH−1G>satisfies Assumtion 5.1. Therefore, they can be solved using Algorithm 11. Steps 3 and 11 of Algorithm 14 perform the assignment zk←min z 1 2z>Hz +q> kz s.t. z ≤z≤z, where qk=(q−G>λif k= 0 q−G>yk−1otherwise, which satisfies Assumtion 5.8, and can therefore be solved using Algorithm 13. Algorithm 16 shows the particularization of Algorithm 14 to the QP problems of the MPC formulations (5.15) and (5.16) under Assumtion 5.18. The control action u(t) is taken as the first melements of ˜z∗. Remark 5.19. One of the key aspects of Algorithm 16 is that the matrix-matrix and matrix-vector operations can be performed sparsely without needing to store all the non-zero elements of the matrices. Instead, it only requires the computation/storage of the repeating elements once, such as matrices Aand B, for G; Q−1,R−1and T−1(if applicable) for H−1;−Qxr,−Rurand −T xrfor q;−Ax(t) for b; and x,x,u,ufor zand z. Remark 5.20. We note that steps 5 and 13 of Algorithm 16 can be solved using Algorithm 13 because the Hessian is diagonal and we are considering box constraints. However, these conditions can be relaxed without greatly compromising the efficiency of the algorithm. In particular, (i) if non-diagonal cost function matrices are considered, and/or (ii) we consider coupled input-state constraints(4.4), then these steps would result in Ndecoupled small-scale, inequality-constrained QP problems, which could be individually solved using interior point or active set methods.
92 Chapter 5. Sparse solvers for model predictive control Algorithm 16: Sparse FISTA solver for standard MPC formulations Require: x(t)∈Rn,xr∈Rn,ur∈Rm,λ∈Rmz,∈R>0 1Compute qwith xrand ur 2Compute bwith x(t) 3t0←1, k←0 4qk←q−G>λ 5zk←solve boxQP(qk;H, z, z) 6Γk← −(Gzk−b) 7∆λk←solve W(Γk) 8yk←λ+ ∆λk 9λk←λ+ ∆λk 10 repeat 11 k←k+ 1 12 qk←q−G>yk−1 13 zk←solve boxQP(qk;H, z, z) 14 Γk← −(Gzk−b) 15 ∆λk←solve W(Γk) 16 λk←∆λk+yk−1 17 tk←1 21 + q1+4t2 k−1 18 yk←λk+tk−1−1 tk (λk−λk−1) 19 until kΓkk∞≤ Output: ˜z∗←zk,˜ λ∗←yk−1 Remark 5.21. The typical choice for the terminal cost function matrix Tis to take it as the solution of the discrete Riccati equation A>TA −T−(A>TB)(R+B>TB)−1(B>T A) + Q=0n×n. This matrix, however, is generally non-diagonal and thus cannot be used in Algorithm 16 unless the last nelements of zkare updated in step 13 by solving a small-scale inequality-constrained QP problem, as discussed inn Remark 5.20. 5.4.2 ADMM-based solver for standard MPC We apply the results of Section 5.3 to the QP problems of the MPC formulations (5.15) and (5.16), whose ingredients are given by (5.17) and (5.18), respectively, under the following assumption. Assumption 5.22. Let Assumtion 4.7 hold and assume that Q∈Sn +,R∈Sm + and T∈Sn +.
5.4. Simple standard MPC formulations 93 Algorithm 17: Sparse ADMM solver for standard MPC formulations Require: x(t)∈Rn,xr∈Rn,ur∈Rm,v0∈Rnz,λ0∈Rnz,ρ∈R>0, p∈R>0,d∈R>0 1Compute qwith xrand ur 2Compute bwith x(t) 3k←0 4repeat 5qk←q+λk−ρvk 6zk+1 ←solve eqQP(qk, b;Hρ, G) 7ˆqk← −ρzk+1 −λk 8vk+1 ←solve boxQP(ˆqk;ρ, z, z) 9λk+1 ←λk+ρ(zk+1 −vk+1) 10 k←k+ 1 11 until kzk−vkk∞≤pand kvk−vk−1k∞≤d Output: ˜z∗←zk, ˜v∗←vk,˜ λ∗←λk It is easy to see that the satisfaction of Assumtion 5.22 implies the satisfaction of Assumtion 5.11 for both QP problems. Therefore, Algorithm 15 can be applied to them. Let us take a closer look at steps 4 and 6 of Algorithm 15. Step 4 of Algorithm 15 requires solving the QP problem (5.13), which due to the block diagonal structure of Hρunder Assumtion 5.22 and to the banded structure of G, satisfies Assumtion 5.3. Therefore, it can be solved using Algorithm 12 as described in Corollary 5.6. Step 6 requires solving the QP problem (5.14), which satisfies Assumtion 5.8 and can therefore be solved using Algorithm 13. Algorithm 17 shows the particularization of Algorithm 15 to the QP problems of the MPC formulations (5.15) and (5.16) under Assumtion 5.22. The control action u(t) is taken as the first melements of ˜v∗. Remark 5.23. We note that the reason why Assumtion 5.22 considers nondiagonal cost function matrices, when Assumtion 5.18 does not, is that step 8 of Algorithm 17 can always be solved using Algorithm 13, whereas in the FISTA-based solver the use of non-diagonal cost function matrices prevented its use. However, Remark 5.20.(ii) still also holds here for the case in which coupled input-state constraints (4.4) are used in place of box constraints. Remark 5.24. Remark 5.19 also applies to Algorithm 17. That is, the matrixmatrix and matrix-vector operations are all performed sparsely thanks to the simple structures of the matrices and only the absolutely necessary variables are stored and computed.
100 Chapter 5. Sparse solvers for model predictive control Algorithm 18: Sparse ADMM-based solver for (5.19) Require: x(t), (xr, ur), v0,λ0,p>0, d>0 1Compute qwith xrand ur 2Compute bwith x(t) 3k←0 4repeat 5ˆqk←q+ (λk ◦−ρvk ◦, P1/2λk f−ρPvk f) 6zk+1 ←solve eqQP(ˆqk, b;ˆ H, G) 7vk+1 ◦←solve boxQP(−ρzk+1 ◦−λk ◦;ρ, v◦, v◦) 8vk+1 f←zk+1 f+ρ−1P−1/2λk f 9if (vk+1 f−c)>P(vk+1 f−c)> r2then 10 vk+1 f←r(vk+1 f−c) q(vk+1 f−c)>P(vk+1 f−c) +c 11 end if 12 λk+1 ◦←λk ◦+ρ(zk+1 ◦−vk+1 ◦) 13 λk+1 f←λk f+ρP1/2(zk+1 f−vk+1 f) 14 k←k+ 1 15 until rp≤pand rd≤d Output: ˜z∗←zk, ˜v∗←vk,˜ λ∗←λk 5.5.2 Computation of admissible ellipsoidal invariant sets Consider a system described by (4.1) subject to Cx ≤c, Du ≤d, (5.35) where C∈Rpx×n,c∈Rpx,D∈Rpu×m,d∈Rpu; and a steady state reference (xr, ur) satisfying the above constraints that we wish the system to converge to. We describe a procedure taken from various articles in the literature (see, for instance, [136] [137, § C.8.1], [132]) for computing a control gain Kand an admissible positive invariant ellipsoidal set E(P, c, r) (see Definition N.7 and Definition N.8) of system (4.1) subject to (5.35) for the control law u(t)=K(x(t)−xr)+ur such that the closed-loop system admissibly converges to xrfor any initial state x(t)∈ E(P, c, r). The procedure is based on solving an optimization problem subject to LMIs. To simplify the procedure, we start by shifting the origin to the reference (xr, ur). That is, let ˆx. =x−xrand ˆu. =u−ur. Then, we compute a gain K and an ellipsoid E(P, 0n, r) (i.e., taking c=xr), such that the closed loop system ˆx(t+ 1) = AKˆx(t), where AK. =A+BK, admissibly converges to the origin for the constraints Cˆx≤ˆc, DKˆx≤ˆ d, (5.36)
5.5. MPC with terminal quadratic constraint 101 where ˆc∈Rpxand ˆ d∈Rpuare given by ˆc=c−Cxr,ˆ d=d−Dur. The implication for invariance given in Definition N.8 reads as ˆx>Pˆx≤r2=⇒ˆx>A> KPAKˆx≤r2, which, applying the S-procedure [23, § B.2], is satisfied if there exists a scalar λ≥0 such that λP −A> KPAK0n 0> nr2(1 −λ)∈Sn+1 ++ , where we note that λmust therefore satisfy λ∈[0,1). This expression can be rewritten as λP 0n 0> nr2(1 −λ)−A> K 0PAK0∈Sn+1 ++ . Applying the Schur complement [23, § A.5.5] leads to: λP 0nA> K 0> nr2(1 −λ)0> n AK0nP−1 ∈S2n+1 ++ , Finally, preand post-multiplying by diag(P−1,1, In) and taking the transformation W. =P−1,Y. =KP−1, leads to the LMI: λW 0nWA>+Y>B> 0> nr2(1 −λ)0> n AW +BY 0nW ∈S2n+1 ++ .(5.37) Next, (5.36) must be satisfied for all ˆx∈ E(P, 0, r). It is well known that max x∈E(P,0n,r)v>x=r√v>P−1v, where v∈Rn[132, § 5.2.2]. Therefore this condition can be imposed by finding P and Ksuch that r2C> jP−1Cj≤ˆc2 j, j ∈Zpx 1, r2(DjK)>P−1(DjK)≤ˆ d2 j, j ∈Zpu 1, where the subindex jindicates the j-th row of the matrix or component of the vector. Taking the above definitions of Wand Y, this can be posed as the LMIs: r2C> jWCj≤ˆc2 j, j ∈Zpx 1,(5.38) "r−2ˆ d2 jDjY Y>D> jW#∈Sn+1 ++ , j ∈Zpu 1.(5.39)
102 Chapter 5. Sparse solvers for model predictive control The procedure is to select a value of r(typically, for convenience, we pick r= 1) and then to solve the following convex optimization problem subject to LMI constraints min W,Y −tr(W) (5.40a) s.t. (5.37),(5.38) and (5.39),(5.40b) for increasing values of λin the range 0 ≤λ < 1 until a feasible solution is found. Finally, Pand Kare recovered from the solution of (5.40). The minimization of −tr(W) is done to maximize the size of the resulting ellipsoid. Remark 5.29. Similar procedures can be used to compute the terminal ingredients of a robust MPC controller that follows formulation (5.19) [16, § V.A]. 5.6 MPC for tracking The solvers presented in the previous sections all consider variations of the standard MPC formulations (4.7). This section, on the other hand, presents a sparse solver, based on the extended ADMM algorithm (Algorithm 3), for a non-standard MPC formulation known as MPC for tracking (MPCT) [138, 139]. This formulation differs from standard MPC formulations (4.7) in the inclusion of a pair of decision variables (xs, us)∈Rn×Rmknown as the artificial reference. Out of the different variants of the MPCT formulation, this dissertation is concerned with the following one, which uses a terminal equality constraint. min x,u, xs,us N−1 X j=0 kxj−xsk2 Q+kuj−usk2 R+kxs−xrk2 T+kus−urk2 S(5.41a) s.t. x0=x(t) (5.41b) xj+1 =Axj+Buj, j ∈ZN−1 0(5.41c) x≤xj≤x, j ∈ZN−1 1(5.41d) u≤uj≤u, j ∈ZN−1 0(5.41e) xs=Axs+Bus(5.41f) x+εx≤xs≤x−εx(5.41g) u+εu≤us≤u−εu(5.41h) xN=xs,(5.41i) where εx∈Rnand εu∈Rmare vectors with arbitrarily small positive components which are added to avoid a possible loss of controllability when the constraints are active at the equilibrium point [139]. We consider the following assumption, which we note does not assume that (xr, ur) is an admissible steady state of system (4.1) subject to (4.3).
5.6. MPC for tracking 103 Assumption 5.30. Let Assumtion 4.2 hold and assume that: (i) Q∈Sn +,R∈Sm +,T∈Sn +and S∈Sm +. (ii) x(t) strictly belongs to the feasibility region of (5.41). That is, there exist x,u,xsand ussuch that (5.41b), (5.41c), (5.41f) and (5.41i) are satisfied and (5.41d), (5.41e), (5.41g) and (5.41h) are strictly satisfied. Remark 5.31. We note that Assumtion 5.30.(ii) implies that x+εx< x −εx and u+εu< u −εu. The cost function (5.41a) penalizes, on one hand, the difference between the predicted states xjand control actions ujwith the artificial reference xsand us, respectively, and on the other, the discrepancy between the artificial reference and the reference (xr, ur) given by the user. The inclusion of the artificial reference provides the MPCT formulation with a series of advantages with respect to other (standard) MPC formulations. First, a common issue of standard MPC formulations with stability guarantees is that the domain of attraction of the controller can become insufficient if the prediction horizon is chosen too small. However, the use of small prediction horizons is desirable in order to help overcome the computational and memory limitations typically imposed by embedded systems. The MPCT formulation provides significantly larger domains of attraction than standard MPC formulations [139], especially for small prediction horizons. Second, it intrinsically deals with references that are not attainable, i.e., that are not a steady state of the system and/or that violate the system constraints [139]. In this case, it will steer the closed-loop system to the admissible steady state (xa, ua)∈Rn×Rmthat minimizes the cost kxa−xrk2 T+kua−urk2 S. Finally, it also guarantees recursive feasibility of the closed-loop system even in the event of a sudden reference change [138]. We now present a sparse solver, originally presented in [10], for the MPCT formulation (5.41) that is based on the EADMM algorithm (Algorithm 3). The MPCT formulation (5.41) can be expressed as a QP problem (5.12). However, the inclusion of the artificial reference (xs, us) leads to a Hessian matrix that is not block diagonal. Therefore, if we were to follow the same procedure presented in Section 5.3, we would not be able to solve step 4 of Algorithm 15 using Algorithm 12 because the resulting Wmatrix would not satisfy Assumtion 5.1. Therefore, we propose to use instead use EADMM (Algorithm 3), which will allow us to recover the structures that emerged in the previous sections, thus resulting in a sparse solver with a very similar iteration complexity to the previous ones.
104 Chapter 5. Sparse solvers for model predictive control 5.6.1 Recasting the MPCT formulation for EADMM Let us start by recasting problem (5.41) by defining the auxiliary variables ˜xi. = xi−xsand ˜ui. =ui−usas follows min ˜ x,˜ u,x, u,xs,us N X i=0 k˜xik2 Q+k˜uik2 R+kxs−xrk2 T+kus−urk2 S(5.42a) s.t. x0=x(t) (5.42b) ˜xi+1 =A˜xi+B˜ui, i ∈ZN−1 0(5.42c) x≤xi≤x, i ∈ZN−1 1(5.42d) u≤ui≤u, i ∈ZN−1 0(5.42e) x+εx≤xN≤x−εx,(5.42f) u+εu≤uN≤u−εu,(5.42g) xs=Axs+Bus(5.42h) ˜xi+xs−xi=0n, i ∈ZN 0(5.42i) ˜ui+us−ui=0m, i ∈ZN 0(5.42j) xN=xs,(5.42k) uN=us,(5.42l) where the decision variables are xs,us,˜ x= (˜x0,...,˜xN), ˜ u= (˜u0,...,˜uN), x= (x0, . . . , xN), and u= (u0, . . . , uN). Equality constraints (5.42i) and (5.42j) impose the congruence of the decision variables with the original problem (5.41). Inequalities (5.41g) and (5.41h) are omitted because they are already imposed by (5.42f) and (5.42g) alongside the inclusion of (5.42k) and (5.42l). Remark 5.32. Note that the summation in the cost function (5.42a) now includes i=N. However, this does not change the solution of the optimization problem due to the inclusion of (5.42k) and (5.42l). We can now obtain a problem of form (2.13) by taking z1= (x0, u0, x1, u1, . . . , xN−1, uN−1, xN, uN),(5.43a) z2= (xs, us),(5.43b) z3= (˜x0,˜u0,˜x1,˜u1,...,˜xN−1,˜uN−1,˜xN,˜uN),(5.43c) which leads to θ1(z1)=0, θ2(z2) = 1 2z> 2diag(T, S)z2−(Txr, Sur)>z2, θ3(z3) = 1 2z> 3diag(Q, R, Q, R, . . . , Q, R)z3,
5.6. MPC for tracking 105 C1= [In0n×m]0 0 −In+m0 0 0...0 0 0 −In+m 0 0 −In+m , C2= 0 In+m . . . In+m In+m , C3= 0. . . 0 In+m0 0 0...0 0 0 In+m 0. . . 0 , b = x(t) 0 . . . 0 0 . Matrices C1,C2and C3contain the equality constraints (5.42b), (5.42i), (5.42j), (5.42k) and (5.42l). Specifically, the first nrows impose constraint (5.42b), the last n+mrows impose the constraints (5.42k) and (5.42l), and the rest of the rows impose the constraints (5.42i) and (5.42j). Set Z1is the set of vectors z1(5.43a) for which the box constraints (5.42d)-(5.42g) are satisfied; set Z2is the set of vectors z2(5.43b) that satisfy the equality constraint (5.42h); and set Z3is the set of vectors z3(5.43c) that satisfy the equality constraints (5.42c). Remark 5.33. Our selection of ziand Cifor i∈Z3 1results in an optimization problem that satisfies Assumption 2.8. Therefore, under a proper selection of ρ, the iterates of the EADMM algorithm will converge to the optimal solution of the MPCT controller. In practice, the parameter ρmay be selected outside the range shown in Theorem 2.9 in order to improve the convergence rate of the algorithm [36]. In this case, the convergence will not be guaranteed and will have to be extensively checked with simulations. 5.6.2 EADMM-based solver for the MPCT formulation We now particularize Algorithm 3 to the optimization problem (5.42). Let us take a closer look at steps 3, 4 and 5 of Algorithm 3. Step 3 of Algorithm 3 minimizes the Lagrangian (2.14) over z1, resulting in the following box-constrained optimization problem: min z1 1 2z> 1H1z1+q> 1z1(5.44a) s.t. z1≤z1≤z1,(5.44b) where H1=ρC> 1C1, q1=ρC> 1C2zk 2+ρC> 1C3zk 3+C> 1λk−ρC> 1b, z1= (−Mn, u, x, . . . , u, x +εx, u +εu), z1= (Mn, u, x, . . . , u, x −εx, u −εu),
106 Chapter 5. Sparse solvers for model predictive control and Mn∈Rn>0 has arbitrarily large components. Due to the structure of C1, we have that H1∈D(N+1)(n+m) ++ . Therefore, problem (5.44) can be solved using Algorithm 13, since it satisfies Assumtion 5.8. Remark 5.34. Remark 5.20.(ii) also applies to problem (5.44). Step 4 of Algorithm 3 minimizes the Lagrangian (2.14) over z2= (xs, us), resulting in the following equality-constrained QP problem: min z2 1 2z> 2H2z2+q> 2z2(5.45a) s.t. G2z2=b2,(5.45b) where H2= diag(T, S) + ρC> 2C2, q2=−(Txr, Sur) + ρC> 2C1zk+1 1+ρC> 2C3zk 3+C> 2λk−ρC> 2b, G2= [(A−In)B], b2=0n. The solution of this problem is given by Proposition 5.5. Note that, since the prediction model (5.41c) is assumed to be controllable, G2is full rank. Moreover, H2∈Sn+m ++ . Therefore, the optimal solution of (5.45) can be obtained by substituting (5.3a) into (5.3b), which leads to the expression z∗ 2=M2q2, where M2=H−1 2G> 2(G2H−1 2G> 2)−1G2H−1 2−H−1 2∈R(n+m)×(n+m)is computed offline and b2does not appear because it is equal to zero. Step 5 of Algorithm 3 minimizes the Lagrangian (2.14) over z3, resulting in the following equality-constrained QP problem: min z3 1 2z> 3H3z3+q> 3z3 s.t. G3z3=b3, where H3= diag(Q,R,Q,R,...,Q,R) + ρC> 3C3, q3=ρC> 3C1zk+1 1+ρC> 3C2zk+1 2+C> 3λk−ρC> 3b, G3= A B −In0··· ··· 0 0 0 0 A B −In··· 0 0 0 0 ··· .........0 0 0 0 ··· 0A B −In0 , b=0Nn,
5.6. MPC for tracking 107 Algorithm 19: Extended ADMM for MPCT Require: z0 2,z0 3,λ0,ρ > 0, > 0 1Update bwith x(t) 2Update q2with xrand ur 3k←0 4repeat 5q1←ρC> 1C2zk 2+ρC> 1C3zk 3+C> 1λk−ρC> 1b 6zk+1 1←solve boxQP(q1;H1, z1, z1) 7q2← −(Txr, Sur) + ρC> 2C1zk+1 1+ρC> 2C3zk 3+C> 2λk−ρC> 2b 8zk+1 2←M2q2 9q3←ρC> 3C1zk+1 1+ρC> 3C2zk+1 2+C> 3λk−ρC> 3b 10 zk+1 3←solve eqQP(q3,0;H3, G3) 11 Γ←3 P i=1 Cizk+1 i−b 12 λk+1 ←λk+ρΓ 13 k←k+ 1 14 until kΓk∞≤,kzk 2−zk−1 2k∞≤,kzk 3−zk−1 3k∞≤ Output: ˜z∗ 1←zk 1, ˜z∗ 2←zk 2, ˜z∗ 3←zk 3,˜ λ∗←λk which, due to the block diagonal structure of H3and the banded structure of G3 satisfies Assumtion 5.3 and can therefore be solved using Algorithm 12. Algorithm 19 shows the particularization of Algorithm 3 to the MPCT formulation (5.42) obtained from the above discussion. The control action u(t) is obtained from the elements u0of the output ˜z∗(5.43a) returned by Algorithm 19. Remark 5.35. It has been shown that the performance of ADMM can be significantly improved by having different values of ρfor different constraints [140, § 5.2], i.e., by considering ρas a diagonal positive definite matrix. In particular, we find that the convergence of Algorithm 19 improves significantly if the equality constraints (5.42b),(5.42k),(5.42l),(5.42j) for i=N, and (5.42i) for i= 0 and i=N, are penalized more than the others. Remark 5.36. We note the operations in Algorithm 19 with matrices Ci,i∈Z3 1, can be performed exclusively using vector-vector operations. This is due to the fact that the matrices C> iCj,i∈Z3 1,j∈Z3 1, are diagonal, with some having a very small number of off-diagonal elements that are easily identified. This is also true for the case in which ρis taken as a diagonal positive definite matrix (Remark 5.35). We maintain the expressions using matrices Ciin Algorithm 19 for simplicity. Remark 5.37. The theoretical upper bound for ρprovided in Theorem 2.9 is easily computable in this case. Indeed, we have that C> 3C3is the identity matrix,
108 Chapter 5. Sparse solvers for model predictive control Figure 5.1: Double reactor and separator system. and therefore its spectral norm is kC> 3C3k= 1. Furthermore, µ3is the minimum eigenvalue of diag(Q, R), which is simple to compute. However, in practice, we find that better results are obtained following Remark 5.35 for larger values of ρthan the one described in Theorem 2.9. Therefore, the convergence of the algorithm will have to be extensively checked with simulations. 5.7 Test Benches This section presents various systems that will be used in the numerical results shown in the following sections. In particular, we present the following systems: a multivariable chemical plant consisting of two reactors and one separator, a ball and plate system, and a series of masses connected by springs. 5.7.1 Chemical plant: double reactor and separator The chemical plant system, depicted in Figure 5.1 and inspired from [141], is a chemical plant in which the two following first-order reaction take place between the reactants A, B and C: A→B,(5.46a) B→C.(5.46b) The reactions occur in two consecutive cylindrical reactors (labeled 1 and 2), which are fed by flows Ff1and Ff2, respectively, of reactants A and B with concentrations cA0and cB0The reactants then flow into a cylindrical separator (labeled 3), where part of the reactants are redirected either to reactor 1 through flow FRor discarded through flow FD. The concentration of the reactants in each cylinder is denoted by cXi, where X∈ {A,B,C}and i∈Z3 1. Each cylinder i∈Z3 1 has a bottom with area Aiand a height of liquid hithat will be determined by the input and output flows of the cylinder, where the output flows Fiare connected as
5.7. Test Benches 109 Table 5.1: Parameters of the double reactor and separator system Parameter Value Units Parameter Value Units A11 m2T0313 K A21 m2kA10−51/s A31 m2kB5·10−61/s ρ1100 kg/m3EA/R-2840 K Cp4 kJ/kg K EB/R-2077 K kv150 kg/m s ∆HA-100 kJ/kg kv250 kg/m s ∆HB-39 kJ/kg kv330 kg/m s αA3.5 - αC0.5 - αB1.1 - αD0.001 - cA01 wt(%) cB00 wt(%) Table 5.2: Operating point of the double reactor and separator system Var. Val. [m] Var. Val. [wt(%)] Var. Val. [wt(%)] Var. Val. [K] h◦ 10.7 c◦ A10.4155 c◦ B10.5480 T◦ 1329 h◦ 20.9 c◦ A20.2581 c◦ B20.6755 T◦ 2333 h◦ 31.33 c◦ A30.2282 c◦ B30.7 T◦ 3323 shown in Figure 5.1. The kinetics of the reactions (5.46) taking place in reactors i∈Z2 1are given by the Arrhenius equations kAi =kAe−EA RTi, kBi =kBe−EB RTi, respectively, where eis Euler’s number, kAand kBare the pre-exponential factor, EAand EBthe activation energies, Rthe universal gas constant, and Tithe temperature of each reactor. The cylinders i∈Z3 1have heating/cooling systems which transfer heats Qi, thus affecting the temperatures and, as a result, the speed of the reactions. The non-linear model (4.5) of the plant is given by: dh1 dt=1 ρA1 (Ff1+FR−F1),(5.47a) dcA1 dt=1 ρA1h1 (Ff1(cA0−cA1) + FR(cAR −cA1)) −kA1cA1,(5.47b) dcB1 dt=1 ρA1h1 (Ff1(cB0−cB1) + FR(cBR −cB1)) (5.47c) −kB1cB1+kA1cA1,(5.47d) dT1 dt=1 ρA1h1 (Ff1(T0−T1) + FR(TR−T1)) + Q1 ρA1H1Cp (5.47e)
116 Chapter 5. Sparse solvers for model predictive control 0 10 20 30 40 50 0 0.2 0.4 0.6 (a) State trajectory. 0 10 20 30 40 50 -0.4 -0.2 0 0.2 0.4 (b) Control action trajectory. 0 10 20 30 40 50 0 0.5 1 1.5 2 (c) Computation times. 0 10 20 30 40 50 0 100 200 300 400 (d) Number of iterations. Figure 5.5: Closed-loop simulation of the ball and plate: Comparison between the ADMMbased solvers. The black dashed/dotted line in figure (a) represents the upper bound on the velocity of the ball. 0 10 20 30 40 50 0 0.1 0.2 0.3 (a) State trajectory. 0 10 20 30 40 50 -1 -0.5 0 0.5 1 (b) Control action trajectory. 0 10 20 30 40 50 0 0.5 1 (c) Computation times. 0 10 20 30 40 50 0 200 400 600 800 (d) Number of iterations. Figure 5.6: Closed-loop simulation of the oscillating masses: Comparison between the ADMM-based solvers.
5.8. Numerical results 117 Formulation laxMPC equMPC ellipMPC MPCT Iterations Average 130.28 130.64 130.2 197.32 Median 127 128 127 197 Maximum 151 151 151 206 Minimum 127 127 127 196 Comp. time Average 1.170 1.109 1.184 1.796 Median 1.089 1.071 1.137 1.762 Maximum 1.905 1.509 1.675 2.137 Minimum 1.081 1.058 1.129 1.749 Table 5.7: Comparison between the ADMM-based solvers for the chemical plant: number of iterations and computation times [ms]. Formulation laxMPC equMPC ellipMPC MPCT Iterations Average 131.58 120.36 114.44 128.34 Median 87 82.5 53 104 Maximum 233 246 252 351 Minimum 83 68 48 91 Comp. time Average 0.866 0.674 0.667 0.854 Median 0.55 0.571 0.479 0.715 Maximum 2.39 1.248 1.574 2.438 Minimum 0.442 0.357 0.26 0.539 Table 5.8: Comparison between the ADMM-based solvers for the ball and plate: number of iterations and computation times [ms]. Formulation laxMPC equMPC ellipMPC MPCT Iterations Average 183.84 265.9 262.52 157.76 Median 182 269 267 107 Maximum 256 352 397 680 Minimum 111 62 63 63 Comp. time Average 0.239 0.286 0.326 0.259 Median 0.201 0.273 0.306 0.17 Maximum 0.667 0.689 0.774 1.031 Minimum 0.165 0.116 0.135 0.113 Table 5.9: Comparison between the ADMM-based solvers for the oscillating masses: number of iterations and computation times [ms].
118 Chapter 5. Sparse solvers for model predictive control We note that the algorithms still require a significant number of iterations even when the optimal solution does not have any active constraints or when close to the reference. These iterations could be reduced by incorporating a warmstart procedure, such as the ones presented in [145] or [146]. For instance, in [10, § 7], we particularize the warmstart procedure from [145, § II] to the EADMM algorithm for solving the MPCT formulation (Algorithm 19). Its particularization to the other solvers would follow similarly. Finally, we note that the results for the ellipMPC formulation are very similar to the ones obtained for the standard MPC formulations (laxMPC and equMPC ), in spite of the addition of the terminal quadratic constraint. This is even the case during the first 7 sample times in Figure 5.6d, during which the terminal quadratic constraint was active. This indicates that the QCQP problem (5.19) solved using ADMM is comparable, in terms of number of iterations and computation times, to solving standard MPC formulations. 5.8.2 Comparison between the ADMMand FISTA-based solvers We now compare the sparse solvers for standard MPC formulations presented in Section 5.4. That is, we compare the ADMM-based solve given in Algorithm 17 with the FISTA-based solve given in Algorithm 16 applied to the standard MPC formulations (5.15) (labeled equMPC) and (5.16) (labeled laxMPC ). We present the results of applying the four solvers to the three systems described in Section 5.7. The values of Q,Rand Nare shown in Table 5.4. We obtain matrix Tby first computing the matrix ˆ Tthat solves the Riccati equation A>ˆ TA −ˆ T−(A>ˆ TB)(R+B>ˆ TB)−1(B>ˆ TA) + Q=06×6, and then taking Tas the diagonal matrix satisfying T(i,i)= n X j=0 ˆ T(i,j). We note that the reason why matrix ˆ Tis not used, as is typically done in MPC, is because the FISTA-based solvers require the cost function matrices to be diagonal (see Assumtion 5.18). We the above Tbecause it is an easy choice. However, in a real setting we must check that this matrix satisfies the Lyapunov equation (5.49). We start the systems at the operating point and set the reference to the admissible steady states provided in Table 5.5. The exit tolerances of the solvers (ADMM and FISTA) are all set to 10−4and the penalty parameter of the ADMM algorithm is selected as ρ= 15 in all cases. We show the number of iterations and computation times at each sample time. The state and input trajectories are nearly indistinguishable to the naked eye when compared to the ones shown in Section 5.8.1. Figure 5.7 shows the comparison between the solvers applied to the chemical plant described in Section 5.7.1. Figure 5.7a shows the computation times and
5.8. Numerical results 119 Formulation laxMPC equMPC FOM ADMM FISTA ADMM FISTA Iterations Average 129.46 1 130.64 1 Median 126 1 128 1 Maximum 151 1 151 1 Minimum 126 1 127 1 Comp. time Average 1.193 0.048 1.222 0.048 Median 1.083 0.031 1.15 0.033 Maximum 2.12 0.382 1.903 0.417 Minimum 1.063 0.029 1.055 0.032 Table 5.10: Comparison between the ADMMand FISTA-based solvers for standard MPC applied to the chemical plant: number of iterations and computation times [ms]. Formulation laxMPC equMPC FOM ADMM FISTA ADMM FISTA Iterations Average 113.78 138.86 120.36 130.1 Median 73 1 82.5 1 Maximum 243 1203 246 1089 Minimum 59 1 68 1 Comp. time Average 0.682 0.789 0.653 0.728 Median 0.453 0.0245 0.433 0.025 Maximum 1.55 6.288 1.587 5.59 Minimum 0.33 0.023 0.357 0.023 Table 5.11: Comparison between the ADMMand FISTA-based solvers for standard MPC applied to the ball and plate: number of iterations and computation times [ms]. Formulation laxMPC equMPC FOM ADMM FISTA ADMM FISTA Iterations Average 193.26 24.24 265.9 26.96 Median 186 1 269 1 Maximum 307 360 352 279 Minimum 102 1 62 1 Comp. time Average 0.249 0.058 0.289 0.06 Median 0.204 0.017 0.275 0.017 Maximum 0.849 0.517 0.762 0.5 Minimum 0.186 0.016 0.119 0.016 Table 5.12: Comparison between the ADMMand FISTA-based solvers for standard MPC applied to the oscillating masses: number of iterations and computation times [ms].
120 Chapter 5. Sparse solvers for model predictive control 0 10 20 30 40 50 0 0.5 1 1.5 2 (a) Computation times. 0 10 20 30 40 50 0 50 100 150 (b) Number of iterations. Figure 5.7: Closed-loop simulation of the chemical plant: Comparison between the ADMMand FISTA-based solvers for the standard MPC formulations. 0 10 20 30 40 50 0 2 4 6 (a) Computation times. 0 10 20 30 40 50 0 500 1000 1500 (b) Number of iterations. Figure 5.8: Closed-loop simulation of the ball and plate system: Comparison between the ADMMand FISTA-based solvers for the standard MPC formulations. 0 10 20 30 40 50 0 0.5 1 (a) Computation times. 0 10 20 30 40 50 0 100 200 300 400 (b) Number of iterations. Figure 5.9: Closed-loop simulation of the oscillating masses system: Comparison between the ADMMand FISTA-based solvers for the standard MPC formulations. Figure 5.7b the number of iterations of each solver. The state and input trajectories are very similar to the ones shown in Figure 5.4. Table 5.10 shows an analysis of the number of iterations and computation times of each solver. Figure 5.8 shows the comparison between the solvers applied to the ball and plate system described in Section 5.7.2. Figure 5.8a shows the computation times and Figure 5.8b the number of iterations of each solver. The state and input trajectories are very similar to the ones shown in Figure 5.5. Table 5.11 shows an analysis of the number of iterations and computation times of each solver. Figure 5.9 shows the comparison between the solvers applied to the oscillating masses system described in Section 5.7.3. Figure 5.9a shows the computation
5.8. Numerical results 121 0 10 20 30 40 50 0 2 4 6 8 (a) Computation times. 0 10 20 30 40 50 0 50 100 150 (b) Number of iterations. Figure 5.10: Closed-loop simulation of the chemical plant: Comparison between solvers applied to the standard MPC formulation with terminal equality constraint. The computation times of qpOASES are shown divided by 10 to be able to appreciate the computation times of the other solvers. times and Figure 5.9b the number of iterations of each solver. The state and input trajectories are very similar to the ones shown in Figure 5.6. Table 5.12 shows an analysis of the number of iterations and computation times of each solver. The results indicate that the FISTA-based solver (Algorithm 16) has a good performance. However, as shown in Figure 5.8 and Table 5.11, the number of iterations can increase significantly if there are active constraints in the optimal solution. The comparison between the number of iterations of the FISTA-based and ADMM-based solvers is somewhat unfair due to the difference between their exit conditions. For instance, in Figure 5.9, the suboptimal solutions obtained with both solvers were at a very similar distance to the optimal solution of the problems at each sample time. In Figure 5.8, however, even though the exit tolerances of both algorithms where set to 10−4, the suboptimal solutions obtained with the FISTA-based solver where up to several orders of magnitude closer to the optimal solutions at each sample time than the ones obtained with the ADMM-based solvers. This may explain the significant increase in the number of iterations of the FISTA-based solvers shown in Figure 5.8 when compared to the mild increase shown in Figure 5.9, even though in both cases there where active constraints in the optimal solutions during the first few sample times. 5.8.3 Standard MPC subject to terminal equality constraint This section compares the ADMM and FISTA-bases solvers for standard MPC formulation (5.15) against other QP solvers from the literature. In particular, we compare Algorithms 17 and 16 with OSQP and qpOASES to solve the QP problem derived from the standard MPC formulation with terminal equality constraint (5.15). The MPC ingredients, solver parameters and references are the same as the ones used in Section 5.8.2. Figure 5.10 shows the comparison between the solvers applied to the chemical plant described in Section 5.7.1. Figure 5.10a shows the computation times and Figure 5.10b the number of iterations of each solver. The state and input trajectories are very similar to the ones shown in Figure 5.4. Table 5.13 shows an
122 Chapter 5. Sparse solvers for model predictive control Solver ADMM FISTA OSQP qpOASES Iterations Average 130.64 1 25 0 Median 128 1 25 0 Maximum 151 1 25 0 Minimum 127 1 25 0 Comp. time Average 1.116 0.04 0.764 72.09 Median 1.067 0.033 0.704 71.99 Maximum 1.841 0.398 2.081 75.93 Minimum 1.054 0.031 0.693 70.79 Table 5.13: Comparison between solvers applied to standard MPC with terminal equality constraint for the chemical plant: number of iterations and computation times [ms]. analysis of the number of iterations and computation times of each solver. The active-set method of the qpOASES solver exits after 0 iterations at every sample time due to the optimal solution of the QP problems never having active constraints. Even so, the initial computations still result in a computation time that is much larger than the ones obtained with the solvers based on first order methods. The OSQP solver always converges in 25 iterations (although if smaller exit tolerances are used the number of iterations grows significantly). Even so, its computation times are similar to the ones obtained by the ADMM solver. The FISTA solver outperforms all others in this case because, once again, it converges in a single iteration if there are no active constraints in the optimal solution. 5.8.4 Standard MPC without terminal constraint This section compares the ADMM and FISTA-bases solvers for standard MPC formulation (5.16) against other QP solvers from the literature. In particular, we compare Algorithms 17 and 16 with OSQP and qpOASES to solve the QP problem derived from the standard MPC formulation without terminal constraint (5.16). The MPC ingredients, solver parameters and references are the same as the ones used in Section 5.8.2. Figure 5.11 shows the comparison between the solvers applied to the ball and plate system described in Section 5.7.2. Figure 5.11a shows the computation times and Figure 5.11b the number of iterations of each solver. The state and input trajectories are very similar to the ones shown in Figure 5.5. Table 5.14 shows an analysis of the number of iterations and computation times of each solver. In this case the qpOASES solver performs a small number of iterations during the first few sample times, since there are active constraints in the optimal solutions of the resulting QP problems. Once again, the iterations drop to 0 as soon as this is no longer the case. The OSQP solver shows a similar trend to FISTA, although with fewer iterations. Nevertheless, the computation times of FISTA are very similar to the ones obtained with OSQP.
5.8. Numerical results 123 0 10 20 30 40 50 0 2 4 6 8 (a) Computation times. 0 10 20 30 40 50 0 500 1000 1500 (b) Number of iterations. Figure 5.11: Closed-loop simulation of the ball and plate: Comparison between solvers applied to the standard MPC formulation without terminal constraint. The computation times of qpOASES are shown divided by 10 to be able to appreciate the computation times of the other solvers. Solver ADMM FISTA OSQP qpOASES Iterations Average 113.78 138.86 73 1.4 Median 73 1 25 0 Maximum 243 1203 525 10 Minimum 59 1 25 0 Comp. time Average 0.618 0.795 1.131 50.75 Median 0.399 0.024 0.41 49.45 Maximum 1.295 6.869 7.288 57.51 Minimum 0.318 0.023 0.399 48.76 Table 5.14: Comparison between solvers applied to standard MPC without terminal constraint for the ball and plate: number of iterations and computation times [ms]. 5.8.5 MPC subject to terminal quadratic constraint This section compares the ADMM-based solver for the MPC formulation subject to terminal quadratic constraints (5.19) to solving the QP problem that arises from considering the same MPC formulation but substituting the terminal ellipsoidal set with the maximal admissible invariant set of the system, which for controllable linear systems is a polyhedron. In particular, we apply both approaches to the oscillating masses system described in Section 5.7.3. We compare the solver against two alternatives: OSQP [35] and FalcOpt [80, 134]. OSQP will be implemented using the maximal admissible invariant set of the system, which is a polyhedron, in place of the ellipsoidal terminal set, thus resulting in a QP problem. The comparison with this solver will allow us to evaluate the computational advantages that can be obtained by using the ellipsoidal set instead of a polyhedral one. FalcOpt is a solver for MPC subject to terminal ellipsoidal constraint that generates the code of the solver following a similar philosophy to Spcies. Therefore, it serves as a comparison of our proposed solver with other state-of-the-art solvers in the literature. However, FalcOpt considers the case of non-linear MPC, does not consider state constraints and the
124 Chapter 5. Sparse solvers for model predictive control 0 10 20 30 40 50 0 0.1 0.2 0.3 0.4 (a) State trajectory. 0 10 20 30 40 50 -1 -0.5 0 0.5 1 (b) Control action trajectory. 0 10 20 30 40 50 0 1 2 3 4 (c) Computation times. 0 10 20 30 40 50 0 100 200 300 400 (d) Number of iterations. Figure 5.12: Closed-loop simulation of the oscillating masses: Comparison between Algorithm 18 for MPC with terminal quadratic constraints and OSQP applied to MPC using the (polyhedral) maximal admissible invariant set. matrix Pof its terminal ellipsoidal constraint must be equal to the terminal cost function T. That is, its terminal ellipsoidal set is given by E(T, xr, r). The prediction horizon Nand the cost function matrices Qand Rare shown in Table 5.4. The ingredients P,cand rdefining the terminal constraint are computed using the procedure described in Section 5.5.2, where the LMI optimization problem is constructed using YALMIP [103] and solved using the SDPT3 solver [144] taking c=xrand r= 1. A feasible solution of problem (5.40) is found for λ= 0.95. The cost function matrix Tis taken as the solution of the Lyapunov function (5.49) for the matrix Kobtained from the procedure used to compute P. We take the penalty parameter of the ADMM algorithm as ρ= 15. We compute the maximal admissible invariant set of the system using the MPT3 toolbox [147], resulting in a polyhedron {x∈Rn:Atx≤bt}described by a matrix Atwith 274 rows, resulting in as many constraints in the QP problem. We solve this QP problem using the OSQP solver. Vector rof the terminal ellipsoidal set of the FalcOpt solver is obtained using a similar procedure to the one described in Section 5.5.2, but forcing P=Tand taking ras a decision variable. We minimize −rto obtain the largest admissible invariant ellipsoid E(T, xr, r) of the system, obtaining r= 3.6552. We use the FalcOpt solver from commit 5ac104c of its GitHub repository [80] with an exit tolerance of 10−3and its “gradients” option set to “casadi”, which uses CasADi [79] (version v3.4.5) to compute the gradients. Figure 5.12 shows the comparison between the two solvers applied to the oscillating masses system described in Section 5.7.3. Figure 5.12c shows the computa-
5.8. Numerical results 125 Solver Algorithm 18 OSQP Iterations Average 262.52 77 Median 267 50 Maximum 397 450 Minimum 63 50 Comp. time Average 0.312 0.747 Median 0.295 0.455 Maximum 0.861 3.944 Minimum 0.183 0.449 Table 5.15: Comparison between solvers applied to MPC with terminal quadratic constraint for the oscillating masses: number of iterations and computation times [ms]. tion times and Figure 5.12d the number of iterations of each solver. Figure 5.12a shows the trajectory of the state p2and Figure 5.12b the trajectory of the control input Fl. Table 5.15 shows an analysis of the number of iterations and computation times of Algorithm 18 and OSQP. The results using FalcOpt are omitted due to them being several orders of magnitude bigger, either because we were not able to fine-tune the solver correctly and/or because its focus on non-linear MPC makes it less efficient at solving the linear MPC problem. The oscillating masses system is described by 6 states and 2 inputs. Even so, its maximal admissible invariant set results in the inclusion of 274 constraints to the QP problem solved by OSQP. This is a reasonable amount of constraints, but the number of can easily become unmanageable, even for average-sized systems, resulting in an inapplicable optimization problem due to the large computation times. In such a case, the computational results shown in Figure 5.12 and Table 5.15 indicate that our proposed solver could result in an applicable controller due to its smaller number of constraints. Furthermore, the computation of the ellipsoidal AIS is computationally attainable even for large-sized systems (for instance, using the approach described in Section 5.5.2), whereas the computation of a polyhedral admissible invariant set can become very challenging. The use of the maximal admissible invariant set has benefits in terms of the domain of attraction of the MPC controller [125, § 4.2], which would be (possibly) reduced with the use of the quadratic constraint. However, the domain of attraction can be enlarged by increasing the prediction horizon, which, due to the sparse nature of the proposed algorithm, may be an acceptable approach. 5.8.6 MPC for tracking This section compares the proposed EADMM-based solver for the MPCT formulation (5.41) (Algorithm 19) against other alternatives. In particular, we compare it against three other alternatives applied to the ball and plate system described in Section 5.7.2. The first two alternatives are the OSQP and qpOASES solvers applied to the
132 Chapter 5. Sparse solvers for model predictive control 012345678 0 200 400 600 800 1000 1200 (a) MPC without terminal constraint. 0123456789 0 100 200 300 400 500 (b) MPC with terminal equality constraint. Figure 5.17: Number of iterations of FISTA applied to the standard MPC formulations using different restart schemes. Closed loop results on the oscillating masses. laxMPC Restart scheme None Alg. 7 Alg. 8 Alg. 10 EfEgE∗ f Avg. Iter. 182.9 107.4 328.3 104.6 235.5 212 338.4 Med. Iter. 114.5 92.5 232 82.5 110 91.5 211 Max. Iter. 472 195 904 228 1166 1101 1151 Min. Iter. 67 39 112 40 38 35 106 equMPC Restart scheme None Alg. 7 Alg. 8 Alg. 10 EfEgE∗ f Avg. Iter. 396.1 506.9 603.5 549 469.2 615.2 651.3 Med. Iter. 400 511 589 553 475 609 633 Max. Iter. 491 816 715 639 622 838 775 Min. Iter. 202 237 513 410 344 440 548 Table 5.18: Analysis of the number of iterations of FISTA with different restart schemes during the sample times shown in Figure 5.17. Application to the oscillating masses.
5.8. Numerical results 133 0 50 100 150 200 250 300 350 400 450 500 10 -4 10 -3 10 -2 10 -1 10 0 (a) Dual norm of the composite gradient mapping. 0 50 100 150 200 250 300 350 400 450 500 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 (b) Distance to the optimum in terms of the objective function value. Figure 5.18: Evolution of the iterates of FISTA for different restart schemes during the first iteration of Figure 5.17a, i.e., applied to laxMPC for the oscillating masses.
134 Chapter 5. Sparse solvers for model predictive control 0 50 100 150 200 250 300 350 400 450 500 10 -4 10 -3 10 -2 10 -1 10 0 (a) Dual norm of the composite gradient mapping. 0 50 100 150 200 250 300 350 400 450 500 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 (b) Distance to the optimum in terms of the objective function value. Figure 5.19: Evolution of the iterates of FISTA for different restart schemes during the first iteration of Figure 5.17b, i.e., applied to equMPC for the oscillating masses.
5.9. Conclusions and future lines of work 135 ures 5.18 and 5.19 show the evolution of the iterates of FISTA for each restart scheme during the first sample time of the closed-loop simulations depicted in Figures 5.17a and 5.17b, respectively. Curiously, the results indicate that the application of restart schemes for the application of FISTA to MPC optimization problems is not always beneficial. In the case of the oscillating masses system, the results show that some of the restart schemes perform better that the non-restarted version, but others provide little to no benefit, with a few moments of significantly worse performance (see Figure 5.17b). The results shown in Table 5.18 show that Algorithms 7 and 10 performed better that the other alternatives, although the extra computational burden required to evaluate their restart conditions may make them less effective in practice than the non-restarted variant. In the case of the ball and plate system, however, the non-restarted version of FISTA performed better that any of the other alternatives. 5.9 Conclusions and future lines of work The numerical results shown in Section 5.8 indicate that the proposed solvers are suitable for their implementation in embedded systems. Their computation times are in the range of milliseconds in all the case studies shown in this dissertation and in the results shown in [12], which were obtained using a Raspberry Pi model 4B. The computation times obtained in the PLC are larger. However, these results were obtained with a less refined version of the solvers and were obtained using a very low resource PLC device. The result show that the proposed solvers, which are applied to relevant MPC formulations, perform well when compared to other state-of-the-art solvers. The application of the restart schemes to solving the MPC optimization problems have, for the most part, not provided benefits in terms of the number of iterations of the algorithm. Additional research is required to determine why this is the case. One possibility is that the matrix Wthat we use leaves little room for improving the non-restarted variant. A future line of work is to study the application of restart schemes to MPC in further detail. The solvers we present are for the most part only applicable in an academic environment. Their application to a real system would require the inclusion of additional ingredients such as a state estimator or a steady state target optimizer. In [8] we included some of these ingredients in the PLC, achieving an offsetfree control of the non-linear system. Additionally, the solvers themselves can be improved further by studying aspects such as: numerical conditioning, the inclusion of soft constraints, certification of the number of iterations, etc. These aspects are all possible future lines of work. The solvers we present here have been included in the SPCIES toolbox [4]. It takes the model of the system, the parameters of the desired MPC formulation and the options of the solver, and automatically generates the library-free code of
136 Chapter 5. Sparse solvers for model predictive control the solver for the target embedded system. Currently, the toolbox only considers the generation of plain C-code and MEX files for Matlab. A future line of work is to maintain the development of this toolbox, including additional programming languages (Python, Julia, etc.) and embedded systems (PLCs, FPGAs, etc.), and adding the above mentioned improvements.
137 Chapter 6 Harmonic based model predictive control for tracking This chapter considers system (4.1) subject to the coupled input-state constraints (4.4). As discussed in Section 4.2, MPC formulations often make use of a terminal admissible invariant set Xt⊆Rn(4.7e) to guarantee the asymptotic stability of the closed-loop system to the desired reference (xr, ur)∈Rn×Rm. However, the use of this terminal set leads to two downsides when the reference to be tracked can change online. The first issue is that the terminal set depends on the value of the reference. If there are a known-before-hand, finite number of references, then a terminal set can be computed offline for each one of them and used online when necessary. Otherwise, the terminal set must be recomputed each time the reference changes, which is typically very computationally demanding. The second issue is that recursive feasibility of the MPC controller can be lost in the event of a reference change, i.e., there may not be a feasible solution of the MPC optimization problem for the current state if the reference has changed since the previous sample time. This second issue is related to the domain of attraction of the MPC controller (Definition 4.6), since feasibility is lost when the current state is no longer within the domain of attraction of the MPC controller due to the reference change. The terminal constraint is the main contributor to this issue when the prediction horizon is not large enough. To see this, note that the predicted state must be able to reach the terminal set within the prediction horizon window, and that systems are typically subject to input constraints. These issues are particularly relevant when dealing with the online implementation of MPC in embedded systems, since their limited computational and memory resources make them unsuitable for large prediction horizons and the on-line computation of admissible invariant sets. One possible solution to avoid having to recompute an admissible positive invariant set is to use a terminal equality constraint, i.e., to employ an MPC
138 Chapter 6. Harmonic based model predictive control for tracking formulation (5.15). However, as discussed in Section 5.4, this formulation severely suffers from small domains of attraction if the prediction horizon is not large enough. Another option is to employ the MPC for tracking formulation (5.41) discussed in Section 5.6. This formulation also uses a terminal equality constraint (5.41i), thus not requiring the computation of a new terminal set each time the reference changes. Moreover, its domain of attraction is significantly larger when compared to other standard MPC formulations, especially for small prediction horizons, although it can be further enlarged if soft constraints are considered [148]. Finally, its recursive feasibility and asymptotic stability are guaranteed even in the event of a reference change. However, as we will show further ahead, the closed-loop performance of this formulation may suffer in certain systems if the prediction horizon is too small. This section presents a novel MPC formulation, originally presented in [13] and [14], which we call harmonic MPC and label by HMPC. The intent behind its development is to obtain an MPC formulation that is suitable for its implementation in embedded systems whilst solving the above issues. In particular, we are interested in an MPC formulation that (i) does not require the computation of a terminal set each time the reference changes, (ii) provides recursive feasibility and asymptotic stability even in the event of a reference change, (iii) has a large domains of attraction and a good closed-loop performance even when using small prediction horizons, and (iv), that the previous three points do not come at the expense of an overly complicated optimization problem, i.e., one that would be more complex to solve online than to simply increase the prediction horizon of the other MPC formulations discussed above. In short, the objective is to develop an MPC formulation with similar properties to the MPCT formulation (5.41) but with an increased domain of attraction and performance for small prediction horizons. In fact, the HMPC formulation can be viewed as an extension of the MPCT formulation (5.41). 6.1 The harmonic MPC formulation The idea behind the HMPC formulation is to substitute the artificial reference (xs, us) of the MPCT formulation (5.41) by the artificial harmonic reference sequences {xhj},{uhj},j∈Z, given by xhj =xe+xssin(w(j−N)) + xccos(w(j−N)),(6.1a) uhj =ue+ussin(w(j−N)) + uccos(w(j−N)),(6.1b) where N∈Z>0is the prediction horizon of the MPC controller and w∈R≥0is the base frequency. The harmonic sequences {xhj}and {uhj}are parameterized by decision variables xe∈Rn,xs∈Rn,xc∈Rn,ue∈Rm,us∈Rm, and uc∈Rm.
6.1. The harmonic MPC formulation 139 To simplify the text, we use the following notation, xH. = (xe, xs, xc)∈Rn×Rn×Rn, uH. = (ue, us, uc)∈Rm×Rm×Rm, ye. =Exe+Fue, ys. =Exs+Fus, yc. =Exc+Fuc,(6.2) where we recall that the matrices E∈Rp×nand F∈Rp×mdefine, along with y and y, the constraint set Y(4.4). For a given prediction horizon Nand base frequency w, the HMPC control law for a given state x(t) and reference (xr, ur) is derived from the following second order cone programming problem labeled by H(x;xr, ur), H(x;xr, ur). = min x,u,xH,uH Jh(x,u,xH,uH;xr, ur) (6.3a) s.t. xj+1 =Axj+Buj, j ∈ZN−1 0(6.3b) y≤Exj+Fuj≤y, j ∈ZN−1 0(6.3c) x0=x(t) (6.3d) xN=xe+xc(6.3e) xe=Axe+Bue(6.3f) xscos(w)−xcsin(w) = Axs+Bus(6.3g) xssin(w) + xccos(w) = Axc+Buc(6.3h) qy2 s(i)+y2 c(i)≤ye(i)−y(i)−εy(i), i ∈Zp 1(6.3i) qy2 s(i)+y2 c(i)≤y(i)−εy(i)−ye(i), i ∈Zp 1,(6.3j) where y∈Rpis a vector with arbitrarily small positive components, and the cost function Jh(x,u,xH,uH;xr, ur) = `h(x,u,xH,uH) + Vh(xH,uH;xr, ur) is composed of two terms: the summation of stage costs `h(·) = N−1 X j=0 kxj−xhjk2 Q+kuj−uhjk2 R,(6.4) and the offset cost function Vh(·) = kxe−xrk2 Te+kue−urk2 Se+kxsk2 Th+kxck2 Th+kusk2 Sh+kuck2 Sh,(6.5) where we consider the following assumption. Assumption 6.1. Let Assumtion 4.2 hold and assume that: (i) y+εy< y −εy. (ii) Q∈Sn ++,R∈Sm ++,T∈Dn ++ and S∈Dm ++.
140 Chapter 6. Harmonic based model predictive control for tracking We denote the optimal value of optimization problem (6.3) for a state x(t) and a given reference (xr, ur) by H∗(x(t); xr, ur)=Jh(x∗,u∗,x∗ H,u∗ H;xr, ur), where x∗, u∗,x∗ H,u∗ Hare the arguments that minimize (6.3). Furthermore, for every j∈Z, we denote by x∗ hj =x∗ e+x∗ ssin(w(j−N)) + x∗ ccos(w(j−N)),(6.6a) u∗ hj =u∗ e+u∗ ssin(w(j−N)) + u∗ ccos(w(j−N)),(6.6b) the harmonic signals parameterized by (x∗ H,u∗ H). At each discrete time instant t, the HMPC control law is given by u(t) = u∗ 0. The use of the harmonic artificial reference is heavily influenced by the extensions of the MPCT formulation to the problem of tracking periodic references [149, 150, 151]. However, in this case, even though the artificial harmonic reference (6.1) is periodic, the reference (xr, ur) to be tracked is a (piecewise) constant set-point. Periodic MPC for tracking formulations are used to track a generic periodic signal with period T∈Z>0, i.e., a reference xr(t), t∈Z, that satisfies xr(t) = xr(t+T)∀t, by employing artificial periodic reference signals {xsj},{usj},j∈Z, with period T. The decision variables xs0to xsT and us0to usT must satisfy the system dynamics (equality constraints) and constraints (inequality constraints). Therefore, the amount of decision variables and constraints grows with T. Let us now focus on the problem of tracking a set-point (xr, ur), and note that this reference can be viewed as a (constant) periodic reference signal whose period is any arbitrary T∈Z>0. Therefore, a periodic MPC for tracking formulation can be used, where the selection of the period Tnow becomes a tuning parameter of the controller. It is rather intuitive that the selection of Twill have an effect on both the performance and domain of attraction of the controller. Its effect on the performance of the controller is not clear, but it seems reasonable to assume that larger periods Twould result in a larger domain of attraction, since we are allowing for more degrees of freedom. However, as discussed above, this would come at the expense of an increase in the number of decision variables and constraints. The idea behind the HMPC formulation is to use an artificial periodic reference signal whose period does not affect the complexity of the optimization problem, and which can therefore be selected to improve the properties of the controller. In particular, the period of the artificial harmonic reference (6.1) is determined by the design parameter w. However, the constraints required to impose the system dynamics and constraints do not depend on the value of w. Indeed, constraints (6.3f)-(6.3h) impose the system dynamics (4.1) and constraints (6.3i)-(6.3j) impose the constraints (4.4), as we prove further ahead (see Corollaries 6.5 and 6.8). Constraint (6.3e) imposes that the predicted state xNreaches the harmonic artificial reference at the end of the prediction horizon, since xhN =xe+xc. Then, noting that the artificial harmonic reference satisfies the system dynamics and constraints, it acts as an admissible invariant set of (4.1) under (4.4). That it, the system would be able to remain in a stable admissible trajectory for all j > N by applying the admissible control actions (6.6b).
6.1. The harmonic MPC formulation 141 Problem (6.3) is a second order cone programming problem due to the inclusion of the constraints (6.3i)-(6.3j), instead of the QP problem typically derived from most linear MPC formulations. However, this class of convex optimization problem is a well studied problem in the literature for which several efficient solvers are available, such as COSMO [152] or ECOS [153]. Remark 6.2. Note that the constraints (6.3b)-(6.3j) do not depend on the reference. Therefore, the feasibility region (Definition 4.5) of the HMPC controller is independent of the reference. As such, feasibility is never lost in the event of reference changes. The following proposition relates constraints (6.3e)-(6.3h) to the satisfaction of the system dynamics. Its proof makes use of the following lemma. Lemma 6.3. Let the elements v`∈Rnvof a sequence {v`}be given by v`=ve+vssin(w`) + vccos(w`),∀`∈Z, where w∈Rand ve∈Rnv,vs∈Rnv, and vc∈Rnv. Then, (i) v`+1 =ve+v+ ssin(w`) + v+ ccos(w`),∀`∈Z, (ii) v+ s(i)2+v+ c(i)2=v2 s(i)+v2 c(i), i ∈Znv 1, where v+ s=vscos(w)−vcsin(w) and v+ c=vssin(w) + vccos(w). Proof: The proof relies on the following well-known trigonometric identities sin(α+β) = sin(α) cos(β) + cos(α) sin(β) cos(α+β) = cos(α) cos(β)−sin(α) sin(β). From these expressions we obtain sin(w(`+ 1)) = sin(w) cos(w`) + cos(w) sin(w`) cos(w(`+ 1)) = cos(w) cos(w`)−sin(w) sin(w`). Therefore, v`+1 =ve+vssin(w(`+ 1)) + vccos(w(`+ 1)) =ve+vs[sin(w) cos(w`)+ cos(w) sin(w`)] +vc[cos(w) cos(w`)−sin(w) sin(w`)] =ve+ [vscos(w)−vcsin(w)] sin(w`)+[vssin(w) + vccos(w)] cos(w`) =ve+v+ ssin(w`) + v+ ccos(w`), which proves claim (i). Denote now Hw. =cos(w)−sin(w) sin(w) cos(w).