Certification of the proximal gradient method under fixed-point arithmetic for box-constrained QP problems

Certification of the proximal gradient method under fixed-point arithmetic for box-constrained QP problems? Pablo Krupaa,∗, Omar Inversob, Mirco Tribastonec, Alberto Bemporadc aSystems Engineering and Automation Department, Universidad de Sevilla, Seville, Spain. bGran Sasso Science Institute (GSSI), L’Aquila, Italy cIMT - School for Advanced Studies, Lucca, Italy Abstract In safety-critical applications that rely on the solution of an optimization problem, the certification of the optimization algorithm is of vital importance. Certification and suboptimality results are available for a wide range of optimization algorithms. However, a typical underlying assumption is that the operations performed by the algorithm are exact, i.e., that there is no numerical error during the mathematical operations, which is hardly a valid assumption in a real hardware implementation. This is particularly true in the case of fixed-point hardware, where computational inaccuracies are not uncommon. This article presents a certification procedure for the proximal gradient method for box-constrained QP problems implemented in fixed-point arithmetic. The procedure provides a method to select the minimal fractional precision required to obtain a certain suboptimality bound, indicating the maximum number of iterations of the optimization method required to obtain it. The procedure makes use of formal verification methods to provide arbitrarily tight bounds on the suboptimality guarantee. We apply the proposed certification procedure on the implementation of a non-trivial model predictive controller on 32-bit fixed-point hardware. Keywords: Convex optimization, Embedded Systems, Predictive control, Fixed-point arithmetic, Gradient method, Certification 1. Introduction Quadratic programming (QP) problems arise in various areas of systems engineering and control, such as model predictive control (MPC), see Rawlings et al. (2017), or reference governors, see Garone et al. (2017), to name a few. Various practical control-related applications, such as the ones listed above, require solving parameter-dependent QP problems at regular intervals on embedded hardware, which poses a challenge due to computational and memory limitations. In recent years there has been a significant advance in this area due to the proposal of efficient QP solvers, some of them for generic QP problems, such as the OSQP solver presented in Stellato et al. (2020), and some tailored to specific problems, such as the solvers proposed in Krupa et al. (2021b); Frison and Diehl (2020), which address MPC optimization problems. In many practical applications of MPC, such as safety-critical systems and space applications, the certification of the maximum number of iterations required by the optimization algorithm and a guarantee of the suboptimality of its provided solution are mandatory for real deployment. Most solvers are based on optimization algorithms with well-known convergence ?This work was supported in part by Grant PDC2021-121120-C21 funded by MCIN/AEI/10.13039/501100011033 and by the “European Union NextGenerationEU/PRTR”, and in part by Grant Margarita Salas (grant number 20122) funded by the Ministerio de Universidades and the European Union (NextGenerationEU). ∗Corresponding author Email addresses: [email protected] (Pablo Krupa), [email protected] (Omar Inverso), [email protected] (Mirco Tribastone), [email protected] (Alberto Bemporad) and suboptimality guarantees. The issue is that these guarantees are typically derived considering ideal conditions, e.g., under the assumption that the mathematical operations performed by the algorithm are error-free; an assumption, however, that is no longer valid when the optimization algorithm is implemented on hardware. This is particularly noticeable on fixed-point hardware, where quantization and round-off errors may lead to significant differences with respect to the “exact” counterpart. The magnitude of this difference depends on the number of fractional bits, which must be selected large enough to provide the required guarantees. In linear-time-invariant (LTI) MPC, the use of explicit MPC (Bemporad, 2019) instead of an iterative solver provides a direct certification of the computation time. However, explicit MPC may require a considerable amount of memory to implement and is only applicable to LTI systems (Bemporad, 2019). In Patrinos et al. (2015), the authors present a dual gradientprojection algorithm for MPC tailored to fixed-point arithmetics. The authors present convergence guarantees and concrete guidelines for selecting the fractional precision to obtain the required suboptimality tolerance. The analysis is done using the notion of the inexact oracle from Devolder et al. (2014), which presents a generic framework for analyzing first-order optimization algorithms in which the oracle provides inexact information. This framework can be used to derive convergence rates when inexact gradient information is available by considering the maximum error when computing the gradient. It has, however, two downsides when applied to fixed-point arithmetic. The first is that it only considers errors in the gradient information, i.e., it considers the other operations performed by the algorithm to be exact, which may not always be the case in fixed-point precision. Second, the convergence results are prearXiv:2303.16786v1 [math.OC] 29 Mar 2023 sented in terms of the average of the iterates of the algorithm (instead of with respect to the current iterate), whose value is not generally available in fixed-point arithmetic, since its computation requires dividing by the current number of iterations. Another approach for analyzing the error propagation is to use affine arithmetic (Fang et al., 2003; Vakili et al., 2013). This framework provides less conservative error bounds than simply considering the worst-case error due to how it handles error propagation in affine operations (addition, substraction and multiplication by a constant), although the error bounds are still conservative in the presence of multiplications between variables. The bounds can be improved by taking a probabilistic approach, as proposed in Fang et al. (2003), at the expense of no longer having a guaranteed certification. This framework was used in Nadales et al. (2022) to certify the minimum number of fractional bits required to satisfy the desired error bound when performing Lipschitz interpolation for data-driven learning-based control. In this case a tight certification could be proposed due to the simplicity of the algorithm, which only required affine operations. However, its application to iterative convex optimization algorithms would provide conservative results, due, precisely, to their iterative nature along with the presence of multiplication operations. A non-conservative approach for analyzing error propagation in fixed-point arithmetic was proposed in Simi´c et al. (2022). The technique allows to check arbitrarily tight error bounds via a bit-vector encoding into integer arithmetics and then into propositional satisfiability, which allows to use mature SAT- based technology. The authors use their procedure to calculate accurate error bounds on a first-order optimization algorithm up to a given number of iterations. However, due to the bit-precise encoding the computational cost of the analysis becomes prohibitive after only a few iterations, even for small optimization problems; this is an issue since first-order methods may require a significant number of iterations to converge. In this article we analyze the implementation of the proximal gradient method (PGM) (Parikh and Boyd, 2013) under fixedpoint arithmetic applied to strongly convex QP problems with box constraints; a choice motivated by its simplicity but practical relevance and by its linear convergence guarantees under exact arithmetic. We provide convergence guarantees in terms of the maximum number of iterations of the algorithm as well as suboptimality guarantees of its output. Our certificate is based on the formal verification procedure presented in Simi´c et al. (2022), which we use to derive arbitrarily tight bounds on the quantities that determine the convergence and suboptimality guarantees. The proposed approach provides a procedure for selecting the minimum fractional precision required to guarantee the given suboptimality and computation-time specifications. Notably, the important difference with respect to Simi´c et al. (2022) is that the verification only needs to analyze a single iteration of the PGM, thus remaining tractable for a wider range of problems. The main features of our approach are: (i) The subptimality guarantees are provided in terms of the output of the algorithm, instead of for the averaged iterates used in Patrinos et al. (2015). (ii) The bounds that determine the suboptimality guarantees are obtained using the formal verification procedure from Simi´c et al. (2022), which allows us to provide a non-conservative error-bound of the gradient computation, in that the exact maximum error committed in the computation of the gradient can be approximated to an arbitrarily large precision. (iii) The certification results are formal guarantees, instead of the probabilistic ones that would be obtained using probabilistic affine arithmetic (Fang et al., 2003) or Monte Carlo analysis (Saracco et al., 2012). We present the application of our approach on a non-trivial MPC problem implemented using 32-bit fixed-point arithmetic. Notation: Given two vectors x, y ∈Rn,x≤(≥)ydenotes componentwise inequalities and hx, yiis their standard inner product. The standard Euclidean norm of a vector x∈Rnis denoted by kxk. =phx, xi. The closed ball of radius r≥0 in Rnis defined as the set Bn r . ={x∈Rn:kxk ≤ r}. The indicator function of a set Cis denoted by IC, i.e., IC(x) = 0 if x∈ C and IC(x) = ∞if x6∈ C. The subdifferential of a function fis denoted by ∂f. We denote by R+the set of strictly positive real numbers. We denote by Rn (p.q)⊂Rnthe set of vectors whose integer and fractional parts are representable using pand qbinary digits, respectively. This notion readily extends to the space of matrices Rn×m, where Rn×m (p.q)⊂Rn×mrepresents the space of matrices whose every element is representable in R(p.q). For any set C ⊆ Rn, we denote C(p.q)={x∈Rn (p.q):x∈ C}. It is obvious that x∈ C(p.q)=⇒x∈ C, but not vice-versa. It is also easy to see that C(p0.q0)⊆ C(p.q)for any p≥p0and q≥q0, and therefore that x∈ C(p0.q0)=⇒x∈ C(p.q). 2. Exact proximal gradient method We consider the class of strongly-convex QP problems min x∈Rn 1 2x>Qx +c>x s.t. ` ≤x≤u, (P) where c, `, u, ∈Rn (p0.q0), with `<u, and Q∈Rn×n (p0.q0), for some finite positive p0and q0. In this article we are interested in finding a suboptimal solution of problem (P), using the proximal gradient method (PGM) (Parikh and Boyd, 2013), for any given realization of the ingredients of problem (P) satisfying Q∈Q⊂Rn×n (p0.q0), c ∈ {c∈Rn (p0.q0):cmin ≤c≤cmax},(1a) `∈ {`∈Rn (p0.q0):`min ≤`≤lmax},(1b) u∈ {u∈Rn (p0.q0):umin ≤u≤umax},(1c) where cmin, cmax, `min, `max, umin, umax ∈Rn (p0.q0),cmin ≤cmax, `max < umin, and Qis a compact set. In particular, we are interested in providing convergence guarantees when implementing the PGM in fixed-point arithmetic. To that end, let us start by recalling the classical “exact” PGM, i.e., its implementation when operating using exact arithmetic. We denote by Pthe set of problems (P) whose ingredients satisfy (1). Let X. ={x∈Rn:`≤x≤u},f:Rn→Rbe given by f(x) = 1 2x>Qx+c>x,x∗∈Rnbe the optimal solution of problem (P) and f∗its optimal value, i.e., f∗=f(x∗). Let L, σ ∈R+be the largest and smallest smoothness and strong convexity parameters of ffor any realization of (P)∈P, i.e., the scalars for which the well-known inequalities f(x)≤f(y) + h∇f(y), x −yi+L 2kx−yk2, f(x)≥f(y) + h∇f(y), x −yi+σ 2kx−yk2, are satisfied ∀x, y ∈Rnfor any given (P)∈P. 2 Algorithm 1: PGM applied to (P)∈P Input: x0∈ X0, 0 < ρ ≤L−1 1For each k= 0,1,2, . . . repeat 2xk+1 ← Tρ(xk) = min{u, max{`, xk−ρ∇f(xk)}} Algorithm 1 shows the PGM applied to a realization of problem (P)∈P, where the min and max operators are taken componentwise and 0 < ρ ≤L−1. It generates a sequence {xk} ∈ X starting at an initial point x0∈ X0, where X0⊆ X is the set of possible initial guesses. Step 2 of the algorithm performs the operator Tρ:X → X given by Tρ(x). = arg min y∈X 1 2ky−(x−ρ∇f(x))k2, which, when particularized to problem (P), is the evaluation of the proximal operator (Parikh and Boyd, 2013) at x−ρ∇f(x). Operator Tρin this setting can also be viewed as the evaluation of the composite gradient mapping (Nesterov, 2013, §2). Furthermore, Algorithm 1 is equivalent in this case to the projected gradient method (Parikh and Boyd, 2013, §4,2). The following theorem recalls the linear convergence of Algorithm 1 when working under “exact” arithmetic. In the following section we will derive a similar result when operating under fixed-point arithmetic. Theorem 1 (Theorem 10.29 from Beck (2017)).Let {xk}be the sequence generated by Algorithm 1 starting at x0∈ X 0ap- plied to (P)∈P. Then, (i) kxk−x∗k2≤(1 −ρσ)kkx0−x∗k2,∀k≥0, (ii) f(xk)−f∗≤1 2ρ(1 −ρσ)kkx0−x∗k2,∀k≥1. 3. Proximal gradient method in fixed-point arithmetic The maximum number of iterations of Algorithm 1 required to guarantee a given suboptimality can be certified using Theorem 1. Indeed, an immediate result of Theorem 1 is that kxk−x∗k2≤, for a given ∈R+, is satisfied for every iteration ksatisfying k≥ log  kx0−x∗k2 log(1 −ρσ).(2) We are now interested in providing an iteration and suboptimality certification when Algorithm 1 is implemented using fixed-point arithmetic. Therefore, let us consider Algorithm 1 when working under fixed-point arithmetic for some predetermined choice of integer and fractional precision (p.q) satisfying p≥p0and q≥q0, where we recall that (p0.q0) is the precision under which the ingredients of (P) are representable. Under this paradigm, we can view the implementation of the PGM algorithm as performing an inexact proximal operator, where the source of inexactness is due to the fixed-point arithmetic and representation of variables. Algorithm 2 shows the implementation of Algorithm 1 under fixed-point arithmetic. It generates a sequence of iterates {ˆxk} ∈ X(p.q)starting from an initial point ˆx0∈ˆ X0 (p.q)⊆ X(p.q). Algorithm 2: PGM applied to (P)∈Punder fixedpoint arithmetic with precision (p.q) Input: ˆx0∈ˆ X0 (p.q), 0 < ρ ≤L−1, ˆ∈R+ (p.q),kmax ∈R+ 1For each k= 0,1,2, . . . repeat 2ˆxk+1 ←min{u, max{`, ˆxk−ˆgρ(ˆxk)}} 3until ˆ d2(ˆxk)<ˆor k≥kmax Step 2 evaluates the operator ˆgρ:X(p.q)→ X(p.q), which is defined as the operator that performs the computation of ρ∇f(·) in the fixed-point paradigm. That is, ˆgρ(ˆxk) returns the result of the evaluation of ρ∇f(ˆxk) when performed under fixed-point arithmetic. Thus, ρ∈R+ (p.q)is an obvious requirement of Algorithm 2. Additionally, the algorithm includes an exit condition given by the satisfaction of the condition ˆ d2(ˆxk)<ˆ, where ˆ d2:X(p.q)→R+ (p.q)∪ {0}is the operator that performs the computation of kˆxk+1 −ˆxkk2in fixed-point arithmetic. This exit condition plays a key role in the suboptimality guarantees provided in this section. The value of ˆgρ(ˆxk) will generally differ from the exact value of the expression ρ∇f(ˆxk) due to the arithmetic errors that occur when using fixed-point arithmetic, thus leading to the source of discrepancy between the sequences generated by Algorithms 1 and 2. The magnitude of this discrepancy, which we formalize in the following definition, will depend on the value of the fractional precision q, with higher values of qobviously leading to smaller errors. Definition 1. Given a choice of q, we denote by Ω∈R+a scalar satisfying kρ∇f(ˆx)−ˆgρ(ˆx)k ≤ Ω,∀ˆx∈ X(p.q),∀(P)∈P. The rest of the operations in Step 2 of Algorithm 2 do not incur in any additional error, as formally stated in the following lemma, so they do not contribute towards the discrepancy between the “exact” and “fixed-point” implementations. Lemma 1. Let v, v, ˆx, ˆy∈Rn (p.q), with v≤v, and consider the set C. ={v∈Rn:v≤v≤v}. Then, the result of the fixedpoint computations min{v, max{v, ˆx+ ˆy}} performed with any precision (ˆp.ˆq)satisfying ˆp≥pand ˆq≥qis the exact Euclidean projection of ˆx+ ˆyonto Cif ˆx+ ˆydoes not result in an overflow. The proof of the lemma is omitted because the claim is a direct result of the fact that the min, max and addition operations do not incur in any error under fixed-point arithmetic as long as there is no overflow in the addition. Corollary 1. Variable ˆxk+1 obtained from Step 2 of Algorithm 2 is the exact Euclidean projection of ˆxk−ˆgρ(ˆxk)onto X, assuming that no overflow occurs during the computations. Therefore, ˆxk∈ X(p.q)⊂ X ,∀k≥0. The reader will note that Lemma 1 is only applicable to Algorithm 2 as long as there is no overflow during its execution, i.e., if the integer precision pis large enough. The certification tool presented in Section 4 can be used to compute the minimum value of prequired to avoid overflow. Thus, we henceforth simply consider that pis chosen so that no overflow occurs during the execution of Algorithm 2. Another useful consequence of Lemma 1 is presented in the following lemma, which states that a scalar Ω satisfying Definition 1 also bounds the error in the computation of ˆxk+1. 3 Lemma 2. Consider Algorithm 2 and let Ωsatisfy Definition 1. Then, kˆxk+1 − Tρ(ˆxk)k ≤ Ω,∀ˆxk∈ X(p.q),∀(P)∈P. Proof. The claim is a direct consequence of Corollary 1 and the fact that the projection operator to non-empty closed convex sets is non-expansive (Ryu and Boyd, 2016, §3.1).  We now present the main result of this section, where we characterize the local linear convergence of Algorithm 2, in terms of the error-bound Ω, when sufficiently far away from the optimal solution. Theorem 2. Let {ˆxk}be the sequence generated by Algorithm 2 applied to a realization of problem (P)∈Pwith starting point ˆx0∈ˆ X0 (p.q)and taking 0< ρ ≤L−1,ρ∈R(p.q). Choose ∈R+ satisfying ρσ > 4Ω, where Ω∈R+is given by Definition 1. Then, as long as kˆxk+1 −x∗k ≥  2, the sequence {ˆxk}satisfies: (i) kˆxk−x∗k2≤1−ρσ 1−4Ω−1k kˆx0−x∗k2,∀k≥0. (ii) f(ˆxk)−f∗≤1−4Ω−1 2ρ1−ρσ 1−4Ω−1k kˆx0−x∗k2,∀k≥1. Proof. See Appendix A. Theorem 2 provides a linear convergence result similar to the one shown in Theorem 1, but where the convergence constant degrades by the factor (1 −4Ω−1). That is, the convergence guarantee worsens as Ω increases and as the desired suboptimality tolerance decreases. Since the theorem only holds as long as kˆxk+1 −x∗k ≥  2, we need to be able to check the satisfaction of this condition during the execution of Algorithm 2. Lemma 3. kˆx− Tρ(ˆx)k ≤ 2kˆx−x∗k,∀ˆx∈ X. Proof. From (Alamo et al., 2019, Property 1.(i)), particularized to our problem formulation and notation, we have that f(Tρ(ˆx)) −f∗≤ρ−1hˆx− Tρ(ˆx),ˆx−x∗i − 1 2ρkˆx− Tρ(ˆx)k2, which along with f(Tρ(ˆx)) −f∗≥0, leads to 1 2kˆx− Tρ(ˆx)k2≤ hˆx− Tρ(ˆx),ˆx−x∗i ≤ kˆx− Tρ(ˆx)k · kˆx−x∗k by making use of the Cauchy-Schwarz inequality.  The previous lemma allows us to guarantee that the condition kˆxk−x∗k ≥  2in Theorem 2 holds for the iterates of Algorithm 2 as long as we can guarantee that kˆxk− Tρ(ˆxk)k ≥ . The following assumption allows us to use the exit condition of Algorithm 2 as a means to guarantee that the condition kˆxk−x∗k ≥  2is satisfied at iteration k. Assumption 1. The exit tolerance ˆof Algorithm 2 satisfies ˆ d2(ˆxk)≥ˆ=⇒ kˆxk− Tρ(ˆxk)k2≥2. In the following section we present a tractable procedure for certifying the satisfaction of this assumption. In practice, we find that one can choose so that the smallest value of ˆsatisfying Assumption 1 is the smallest positive representable number in precision q, i.e., 2−q. In this case the exit condition of Algorithm 2 becomes ˆ d2(ˆxk) = 0. Under Assumption 1, the convergence guarantee provided in Theorem 2 holds as long as the exit condition ˆ d2(ˆxk)<ˆis not satisfied. If it is satisfied at some iteration k, then the convergence guarantee provided by Theorem 2 is only guaranteed to hold until iteration k−1. The following theorem provides the suboptimality guarantees of Algorithm 2 when the exit condition is satisfied at some iteration k. The result makes use of the bounds provided in the following definition. The following section will provide computationally tractable procedures for computing arbitrarily tight values of said bounds. Definition 2. Consider Algorithm 2. For a given choice of fractional precision qand tolerance ˆ, we denote by δ, ω, Θ∈R+ the scalars satisfying kˆxk− Tρ(ˆxk)k ≤ δ, kˆxk+1 − Tρ(ˆxk)k ≤ ω, kˆxk−ˆxk+1k ≤ Θ for all ˆxk∈ X(p.q)satisfying ˆ d2(ˆxk)<ˆ,∀(P)∈P. Theorem 3. Let {ˆxk}be the sequence generated by Algorithm 2 applied to a realization of problem (P)∈Pwith starting point ˆx0∈ˆ X0 (p.q)and taking 0< ρ ≤L−1,ρ∈R(p.q). Denote T. =σ−1(ρ−1+L). Then, if ˆ d2(ˆxk)<ˆ, (i) kˆxk+1 −x∗k ≤ ω+δT, (ii) f(ˆxk+1)−f∗≤ρ−1(Θ + Ω)(ω+δT) + 1 2Θ2. Proof. From (Nesterov, 2013, Lemma 3), particularized to our notation, we have that kˆxk− Tρ(ˆxk)k ≥ T−1kTρ(ˆxk)−x∗k, which leads to kˆxk− Tρ(ˆxk)k ≤ δ=⇒ kTρ(ˆxk)−x∗k ≤ T δ. Claim (i) follows from adding the previous inequality with kˆxk+1 − Tρ(ˆxk)k ≤ ωand applying the triangle inequality. By the same procedure, we also derive kˆxk−x∗k ≤ δ(T+ 1). Claim (ii) then follows from particularizing Lemma 5.(i) to y=x∗, using the Cauchy-Schwarz inequality and then taking the previous inequalities along with the inequalities presented in Definitions 1 and 2.  The following corollary gathers the guarantees that are obtained from Algorithm 2 in terms of the error-bounds and tolerances presented throughout this section. Corollary 2 (Suboptimality guarantee of Algorithm 2).Let kmax ≥log 2 4D/log(C), where C. =1−ρσ 1−4Ω−1and Dsatisfies D≥maxˆx∈ˆ X0 (p.q)kˆx−x∗k2. The following hold: (i) If ˆ d2(ˆxk)≥ˆfor all k={0,...kmax}, then kˆxkmax −x∗k ≤  2and f(ˆxkmax )−f∗≤(2−4Ω) 8ρ. (ii) If ˆ d2(ˆxk)<ˆthen kˆxk+1 −x∗k ≤ ω+δT and f(ˆxk+1)−f∗≤ρ−1(Θ + Ω)(ω+δT) + 1 2Θ2. Remark 1. From Lemma 2 we have that the bound ωfrom Definition 2 can be substituted by Ωfrom Definition 1. We find that for high fractional precision there can be an insignificant difference between the two quantities, thus not meriting the additional computation time required to compute ω. 4 Algorithm 3: Example of input algorithm for formal verification Parameters: (p.q), b∈R(p.q), ˆa∈R+,ξ∈R,χ∈R+ Non-deterministic: a∈ A ={a∈Rm:||a||∞≤ˆa} 1r=ha, ai 2µ=br 3if |err(r)|> χ then FAIL 4if exact(µ)< ξ then FAIL 4. Obtaining error-bounds for Algorithm 1 This section presents procedures for obtaining arbitrarily tight values of the bounds Ω, δ,ωand Θ introduced in Definitions 1 and 2 as well as a procedure for checking the satisfaction of Assumption 1. The procedures consist on the application of the formal verification procedure presented in Simi´c et al. (2022) to check the conditions provided in the definitions and assumption. We first introduce the formal analysis technique for fixed-point arithmetic, and then show how it fits within our certification process. 4.1. Formal verification for fixed-point arithmetic Simi´c et al. (2022) consider the problem of estimating the numerical accuracy of algorithms in fixed-point arithmetic with variables of arbitrary precision and possibly non-deterministic values. The idea is to re-compute in a greater precision the result of each fixed-point operation, so that the numerical error can be estimated based on the difference between the two values; at the same time, the different errors are in turn accounted for and propagated through the re-computations. When sufficient precision is used to store the re-computed values, this yields an accurate error tracking for each variable at any point of the computation. The techniques relies on a bit-precise encoding to transform the sequences of operations under analysis into operations in integer arithmetic over vectors of bits; these are in turn encoded as a SAT formula that is satisfiable if and only if the algorithm under analysis exceeds a given bound on the numerical error. The technique is quite accurate in that it allows to formally verify arbitrarily tight bounds on the numerical error up to a given number of iterations. On the other hand, as observed by the authors, the required program unfolding pass along with the bit-vector encoding ends up introducing considerable overhead; the analysis can become quickly intractable, even for small problems and a few iterations. We now provide an illustrative example to introduce the concepts and notation relevant to this article. Consider Algorithm 3 working under a given fixed-point precision (p.q), where ris a fixed-point variable, exact{r}is its “exact” counterpart and err{r}is its error, i.e., exact{r}=r+err{r}. The same applies for the other variables, such as µ, whose exact value will generally differ from its value computed in fixedpoint arithmetic due to the multiplication operations. The error of ris propagated when computing the error of µ, i.e., exact{µ}will contain the value of ha, ai. We can perform assessments on the fixed-point variables, their errors and their exact values (as shown in Steps 3 and 4) for all possible values of the non-deterministic inputs (variable ain this case). If all the assessments are satisfied for all possible values of the nondeterministic variables then the procedure will return a PASS. Otherwise, it will return a FAIL. Algorithm 4: Algorithm for asserting Assumption 1 Parameters: (p.q), Q, u, l, cmin, cmax, ˆ,2 Non-deterministic: ˆx0∈ X(p.q), (P)∈P 1ˆxk+1 ←min{u, max{`, ˆxk−ˆgρ(ˆxk)}} 2d← hˆxk−ˆxk+1,ˆxk−ˆxk+1i 3if d≥ˆand exact{d}< 2then FAIL Algorithm 5: Algorithm for deriving Ω Parameters: (p.q), Q, u, l, cmin, cmax, Ω2 Non-deterministic: ˆx0∈ X(p.q), (P)∈P 1ˆg←ρ(Qˆx0+c) 2e←err{ˆg} 3v← he, ei 4if exact{v}>Ω2then FAIL Algorithm 6: Algorithm for deriving Def. 2 bounds Select: b∈ {δ2, ω2,Θ2} Parameters: (p.q), Q, u, l, cmin, cmax, ˆ Non-deterministic: ˆx0∈ X(p.q), (P)∈P 1ˆxk+1 ←min{u, max{`, ˆxk−ˆgρ(ˆxk)}} 2if hˆxk−ˆxk+1,ˆxk−ˆxk+1i<ˆthen 3switch bdo 4case δ2do s←ˆxk−ˆxk+1 5case ω2do s←err{ˆxk+1} 6case Θ2do 7err{ˆxk+1} ← 0 8s←ˆxk−ˆxk+1 9end case 10 end switch 11 end if 12 v← hs, si 13 if exact{v}> b then FAIL Example 1. We run the verification procedure presented in Simi´c et al. (2022) on Algorithm 3 with p= 8,q= 8,m= 20, ˆa= 0.125,b= 1.5,ξ= 0, and χ= 0.069580078125. We obtain a PASS, that is, there is no value of ainside the box defined by ˆafor which the assertions stated in Steps 3 and 4 of Algorithm 3 are violated. The result of this example highlights one of the main benefit of this procedure, which is the bound χ= 0.069580078125. Variable ris the inner product of aby itself. As stated in Patrinos et al. (2015), the standard theoretical bound for the maximum error committed by an inner product is given by err(ha, ai)≤2−qm, which is equal to 0.078125 for the values of qand min this example. However, the formal verification procedure has found that this theoretical bound can be improved to err(ha, ai)=0.069580078125, which is a 12.28% improvement. That is, the procedure may lead to tighter bounds on the errors committed by the fixed-point algorithm than the ones obtained by simply evaluating its execution under the theoretical worst-case scenario. Interestingly, considering the dimension of a(m= 20), the value of q, and the size of the box, thanks to our specific formulation of the problem at hand we were able to obtain this verification verdict in a very small amount of time and with limited hardware resources (5ms using a standard machine with an Intel i5 processor). 5 4.2. Verification procedure for Algorithm 2 The application of the verification tool presented in the previous subsection to Algorithm 4 certifies if the given ˆand 2 satisfy Assumption 1. Its application to Algorithms 5 certifies if the provided value of Ω satisfies the condition presented in Definition 1. A PASS indicates that the given Ω2satisfies the condition kTρ(ˆxk)−ˆxk+1k2≤Ω2,∀ˆxk∈ X(p.q),∀(P)∈P, whereas a FAIL indicates that there is at least one combination of ˆxk∈ X(p.q)and (P)∈Pfor which the condition is not satisfied. Note that the condition is asserted with respect to Ω2, since the verification tool does not allow the use of the squareroot operation. Similarly, for a given value of ˆ, Algorithm 6 is used to certify if the bounds δ,ωor Θ satisfy the conditions presented in Definition 2. Arbitrarily tight values of Ω, δ,ωor Θ can be obtained by applying the bisection method to Algorithms 5 and 6. Note that Algorithms 4, 5 and 6 only execute a single iteration of Algorithm 2. Thus, the proposed verification procedure remains tractable for moderately-sized problems, in contrast with the approach taken in Simi´c et al. (2022). Remark 2. The verification tool from Simi´c et al. (2022) can be configured to return a FAIL in the event of a numerical overflow. Thus, we can verify that no overflow occurst in Algorithm 2 for a given choice of the integer precision pif the tool does not fail due to an overflow when applied to Algorithm 4. 5. Numerical case study We apply the verification procedures presented in the previous section to certify the fixed-point implementation of the PGM to solve the optimization problem of a linear MPC controller for a discrete-time, time-invariant system given by a state-space model ˜x(t) = A˜x(t) + B˜u(t), where ˜x(t)∈Rnxand ˜u(t)∈Rnuare the state and control input at sample time t. In particular, we consider the system of three masses connected by springs presented in (Krupa et al., 2021a, §3), where we take the mass all three objects equal to 1kg and the spring constants as 1N/m. The 6-dimensional system state is given by the position and velocity of each of the three objects, while the control input is given by the two external forces applied to the outer objects. We take the following MPC formulation: min Np−1 X i=0 k˜xi−˜xrk2 Wx+k˜ui−˜urk2 Wu+k˜xNp−˜xrk2 P(3a) s.t.˜x0= ˜x(t) (3b) ˜xi+1 =A˜xi+B˜ui, i ∈ZNp−1 0(3c) ˜ui= ˜uNc−1, i ∈ZNp−1 Nc(3d) ˜u−≤˜ui≤˜u+, i ∈ZNc−1 0,(3e) where Nc∈R+is the control horizon; Np≥Ncis the prediction horizon; Wx, P ∈Rnx×nxand Wu∈Rnu×nuare positive definite; (xr, ur) are the state and input references; and ˜u−,˜u+∈Rnusatisfying ˜u−≤˜u+define the bounds on the control input. We take Nc= 2, Np= 5, ˜u+= (0.5,0.5), ˜u−=−˜u+,Wx= 0.5Inx,Wu= 0.25Inuand Pas the solution of the associated discrete algebraic Riccati equation. Problem (3) can be transformed into (P) by eliminating the states and rewriting it in condensed form, see e.g., Richter (2012); Jerez et al. (2011), leading to a (nuNc)-dimensional QP problem. In this case, ingredients Q,`and uof problem (P) are fixed, whereas the value of cwill depend on the value of the reference (xr, ur) as well as the current state ˜x(t). We compute cmin and cmax by assuming that the position of the objects belong to the interval [−0.5,0.5]m and the velocities to [−1,1]m/s. A non-deterministic Qwould be taken if we allowed the possibility of changing the weight Wuonline or if we considered a time-varying model of the system. We now certify the PGM applied to the resulting condensed MPC problem when implemented in fixed-point arithmetic on a 32-bit device, where we take p= 10 for the integer precision and q= 21 for the fractional precision (the remaining bit is used for storing the sign). We store the matrices of the QP problem in the selected precision. The resulting problem has L= 4.9645 and σ= 0.3532. We take ρas the largest number representable in R(p.q)that satisfies ρ≤1/L. In our formal verification procedure we used the prototype tool of Simi´c et al. (2022) for generating the bit-vector encoding, CBMC 5.4 (Clarke et al., 2004) for generating the SAT formula from the bit-vector encoding and MiniSat (E´en and S¨orensson, 2004) to check for the satisfiability of the SAT formula. All computations are performed on an Intel i5 processor running at 1.6GHz. We start by computing the value of Ω by selecting an initial value of Ω2and then applying the bisection method on the verification of Algorithm 5 with an exit tolerance of 10−14, i.e., until the difference between the largest and smallest values of Ω2resulting in a FAIL and a PASS, respectively, is smaller than 10−14. Table 1 shows the value of Ω obtained from this procedure, along with the selected exit tolerance, initial guess of Ω2, number of tests resulting in a PASS, number of tests resulting in a FAIL, average computation times of calls resulting in a PASS or FAIL, and total computation time of the bisection method. We take ˆ= 2−q, which is the smallest value it can take, and then find the largest value of satisfying Assumption 4 by applying the bisection method on Algorithm 4. The results are presented in Table 1, where we note that the value of satisfies  > 4Ω/(ρσ)=9.6195 ·10−5. Therefore, we can use the exit condition ˆ d2(ˆxk) = 0. Finally, we obtain the bounds δand Θ following the same bisection procedure used to compute Ω. The results are also presented in Table 1. We take ω= Ω, as stated in Remark 1. Plugging the results into Corollary 2.(ii), we obtain the following: kmax = 250, if ˆ d2(ˆxk) = 0 then kˆxk+1 −x∗k ≤ 0.0389 and f(ˆxk+1)−f∗≤2.7165 ·10−4, are the best suboptimality bounds that can be guaranteed, since the ones from Corollary 2.(i) are smaller. The value of kmax required to obtain the same /2-suboptimality under exact arithmetic is 217, c.f., (2). 6. Conclusions This article has presented a procedure for certifying the implementation of the PGM under fixed-point arithmetic when applied to strongly-convex box-constrained QP problems. We have proven that the PGM maintains a linear convergence guarantee when sufficiently far away from the optimal solution, indicated by the choice of , whose value can be reduced up to a maximum bound given by the fixed-point error-bound. We have then presented a procedure based on recent formal verification tools to obtain a arbitrarily tight values of this errorbound and the other bounds that characterize the suboptimality of the output of the PGM. Finally, we have shown that the computation times of the proposed verification procedures are tractable for a non-trivial MPC example. 6 Bound bValue b2Tol. # P/F Av. PASS time [s] Av. FAIL time [s] Total time [s] Ω 1.711 ·10−62−2q10−14 4/8 218.1 982.8 8738.8 6.8949 ·10−4(4Ω/(ρσ))210−911/4 73.8 75.2 1117.4 δ1.383 ·10−32−q10−910/4 84.7 90.8 1215.0 Θ 1.381 ·10−3δ210−93/8 12.1 28.3 266.6 Table 1: Bounds obtained from the verification procedures for the three-mass-spring case study. Appendix A. Proofs and auxiliary lemmas We start by providing two lemmas whose results are used in the proofs of Theorems 2 and 3. Lemma 4. Consider Algorithm 2. For any α∈R+, αˆxk−ˆxk+1 −ˆgρ(ˆxk)∈∂IX(ˆxk+1),∀k≥0. Proof. The claim follows from ˆxk+1 being the Euclidean projection of ˆxk−ˆgρ(ˆxk) onto X(see Corollary 1) along with the optimality condition of the projection operator (Bertsekas, 2009, Prop. 5.4.7) and the equivalence between the subdifferential of the indicator function of a non-empty convex set and its normal cone (Bertsekas, 2009, Example 5.4.1).  The following lemma particularizes (Alamo et al., 2019, Property 1) to the fixed-point PGM paradigm. Lemma 5. Consider Algorithm 2 and let vk∈ Bn Ωbe the vectors that satisfy ρ∇f(ˆxk) = ˆgρ(ˆxk)+vkfor every k≥0. Denote sk. = ˆxk−ˆxk+1. Then, f(ˆxk+1)−f(y)≤ρ−1hsk+vk,ˆxk+1 −yi+1 2ρkskk2(i) =ρ−1hsk,ˆxk−yi − 1 2ρkskk2+ρ−1hvk,ˆxk+1 −yi(ii) =1 2ρkˆxk−yk2−1 2ρkˆxk+1 −yk2+ρ−1hvk,ˆxk+1 −yi(iii) Proof. From Lemma 4 we have that ρ−1(ˆxk−ˆxk+1 −ˆgρ(ˆxk)) ∈ ∂IX(ˆxk+1). Therefore, from the definition of the subdifferential (Parikh and Boyd, 2013, §2.3), we have that IX(y)≥ IX(ˆxk+1) + ρ−1hˆxk−ˆxk+1 −ˆgρ(ˆxk), y −ˆxk+1i, where taking y∈ X and recalling that ˆxk+1 ∈ X (see Corollary 1), leads to 0≥ρ−1hˆxk−ˆxk+1 −ˆgρ(ˆxk), y −ˆxk+1i.(A.4) From the convexity of fwe have that f(y)≥f(ˆxk) + h∇f(ˆxk), y −ˆxki.(A.5) Additionally, from the L-smoothness of fand since ρ≤L−1, we have that f(ˆxk)≥f(ˆxk+1)− h∇f(ˆxk),ˆxk+1 −ˆxki − 1 2ρkskk2.(A.6) Claim (i) follows from adding (A.4), (A.5) and (A.6) along with the definition of vk. Claims (ii) and (iii) then follow from simple algebraic manipulations; c.f. Property 1.(i) in Alamo et al. (2019).  We now present the proof of Theorems 2, which closely follows the proofs of (Beck, 2017, Theorem 10.16 and Theorem 10.29), although various modifications have to be made to extend the results to the fixed-point arithmetic paradigm. Proof of Theorem 2. Let sk. = ˆxk+1 −ˆxkand γk. = ˆxk−x∗. Consider the function ψ:X → Rgiven by ψ(y) = f(ˆxk) + h∇f(ˆxk), y −ˆxki+IX(y) + 1 2ρky−ˆxkk2. Since ψis an ρ−1-strongly convex function, it follows from (Beck, 2017, Theorem 5.24) that ψ(y)−ψ(ˆxk+1)≥ hµ, y −ˆxk+1i+1 2ρky−ˆxk+1k2,(A.7) ∀ˆxk+1 ∈ X(p.q),∀µ∈∂ψ(ˆxk+1). Since ρ≤L−1, we have that f(·) satisfies the well-known descent lemma (Beck, 2017, Lemma 5.7) f(ˆxk+1)≤f(ˆxk) + h∇f(ˆxk),ˆxk+1 −ˆxki+1 2ρkˆxk+1 −ˆxkk2, ∀ˆxk+1 ∈ X(p.q),∀y∈ X, which along with the definition of ψ and noting that IX(ˆxk+1) = 0 by virtue of Corollary 1, leads to ψ(ˆxk+1) = f(ˆxk) + h∇f(ˆxk), ski+L 2kskk2≥f(ˆxk+1), ∀ˆxk+1 ∈ X(p.q). Thus, we can rewrite (A.7) as ψ(y)−f(ˆxk+1)≥ hµ, y −ˆxk+1i+1 2ρky−ˆxk+1k2,(A.8) ∀ˆxk+1 ∈ X(p.q),∀µ∈∂ψ(ˆxk+1). The subdifferential of ψevaluated at ˆxk+1 is given by ∂ψ(ˆxk+1) = ∇f(ˆxk) + ρ−1(ˆxk+1 −ˆxk) + ∂IX(ˆxk+1).(A.9) From Definition 1 we have that for each k≥0 there exists a vector vk∈ Bn Ωsatisfying ρ∇f(ˆxk) = ˆgρ(ˆxk) + vk. Therefore, we can rewrite (A.9) as ∂ψ(ˆxk+1) = ρ−1(vk+ ˆgρ(ˆxk) + ˆxk+1 −ˆxk) + ∂IX(ˆxk+1). From Lemma 4 we have 0 ∈ρ−1(ˆgρ(ˆxk)+ˆxk+1−ˆxk)+∂IX(ˆxk+1), thus ρ−1vk∈∂ψ(ˆxk+1). This allows us to rewrite (A.8) as ψ(y)−f(ˆxk+1)≥ hLvk, y −ˆxk+1i+1 2ρky−ˆxk+1k2, ∀ˆxk+1 ∈ X(p.q), for some vk∈ Bn Ω. Undoing the expression of ψ(y) and particularizing to y=x∗∈ X leads to f(ˆx∗)−f(ˆxk+1)≥1 2ρkγk+1k2−1 2ρkγkk2+f(ˆx∗)−f(ˆxk) − h∇f(ˆxk),ˆx∗−ˆxki+ρ−1hvk,ˆx∗−ˆxk+1i. 7 Since fis a σ-strongly convex function, we have that (Beck, 2017, Theorem 5.24.(ii)) f(ˆx∗)−f(ˆxk)− h∇f(ˆxk),ˆx∗−ˆxki ≥ σ 2kˆx∗−ˆxkk2. Thus, f(ˆx∗)−f(ˆxk+1)≥1 2ρkγk+1k2−ρ−1−σ 2kγkk2−ρ−1hvk, γk+1i. (A.10) By definition of ˆx∗, we have that f(ˆx∗)−f(ˆxk+1)≤0. Therefore, the right hand side of (A.10) must also be less or equal to 0, which leads to 1 2kˆxk+1 −ˆx∗k2≤1−ρσ 2kˆxk−ˆx∗k2+hvk,ˆxk+1 −ˆx∗i (∗) ≤1−ρσ 2kˆxk−ˆx∗k2+kvkk · kˆxk+1 −ˆx∗k (∗∗) ≤1−ρσ 2kˆxk−ˆx∗k2+ 2Ω−1kˆxk+1 −ˆx∗k2, where in (∗) we are making use of the Cauchy-Schwarz inequality and in (∗∗) of the fact that kvkk ≤ Ω and kˆxk+1 −ˆx∗k> 2. Since by construction ρσ ≤1, the assumption 4Ω−1< ρσ implies 2Ω−1<1 2. Thus, we derive kˆxk+1 −ˆx∗k2≤1−ρσ 1−4Ω−1kˆxk−ˆx∗k2, which leads to claim (i) if applied recursively, where we note that the assumption 4Ω−1< ρσ guarantees that the sequence is convergent. We now prove claim (ii) by rearranging (A.10) and proceeding as follows: f(ˆxk+1)−f(ˆx∗)≤ρ−1−σ 2kγkk2−1 2ρkγk+1k2+ρ−1hvk, γk+1i ≤ρ−1−σ 2kγkk2−1 2ρkγk+1k2+ 2ρ−1Ω−1kγk+1k2 ≤ρ−1−σ 2kγkk2+ρ−12Ω−1−1 2kγk+1k2 (∗) ≤1−ρσ 2ρkγkk2=1−4Ω−1 2ρ1−ρσ 1−4Ω−1kγkk2 (∗∗) ≤1−4Ω−1 2ρ1−ρσ 1−4Ω−1k+1 kˆx0−ˆx∗k2, ∀ˆxk+1 ∈ X(p.q),∀k≥0, where (∗) holds since 2Ω−1<1 2and (∗∗) follows from claim (i). References Alamo, T., Krupa, P., Limon, D., 2019. Restart FISTA with global linear convergence, in: 2019 18th European Control Conference (ECC), IEEE. pp. 1969–1974. Beck, A., 2017. First-Order Methods in Optimization. Society for Industrial and Applied Mathematics. Bemporad, A., 2019. Explicit Model Predictive Control. Springer London. Bertsekas, D.P., 2009. Convex Optimization Theory. Athena Scientific. Clarke, E., Kroening, D., Lerda, F., 2004. A tool for checking ANSI- c programs, in: Jensen, K., Podelski, A. (Eds.), Tools and Algorithms for the Construction and Analysis of Systems. Springer Berlin Heidelberg. volume 2988, pp. 168–176. Series Title: Lecture Notes in Computer Science. Devolder, O., Glineur, F., Nesterov, Y., 2014. First-order methods of smooth convex optimization with inexact oracle. Mathematical Programming 146, 37–75. E´en, N., S¨orensson, N., 2004. An extensible SAT-solver. Lecture notes in computer science 2919, 502–518. Fang, C., Rutenbar, R., Tsuhan Chen, 2003. Fast, accurate static analysis for fixed-point finite-precision effects in DSP designs, in: International Conference on Computer Aided Design, IEEE. pp. 275–282. Frison, G., Diehl, M., 2020. HPIPM: a high-performance quadratic programming framework for model predictive control. IFAC- PapersOnLine 53, 6563–6569. Garone, E., Di Cairano, S., Kolmanovsky, I., 2017. Reference and command governors for systems with constraints: A survey on theory and applications. Automatica 75, 306–328. Jerez, J.L., Kerrigan, E.C., Constantinides, G.A., 2011. A condensed and sparse QP formulation for predictive control, in: IEEE Conference on Decision and Control and European Control Conference, IEEE. pp. 5217–5222. Krupa, P., Jaouani, R., Limon, D., Alamo, T., 2021a. A sparse ADMM-based solver for linear MPC subject to terminal quadratic constraint. arXiv:2105.08419 . Krupa, P., Limon, D., Alamo, T., 2021b. Implementation of model predictive control in programmable logic controllers. IEEE Transactions on Control Systems Technology 29, 1117–1130. Nadales, J.M., Manzano, J.M., Barriga, A., Limon, D., 2022. Efficient FPGA parallelization of lipschitz interpolation for real-time decision-making. IEEE Transactions on Control Systems Technology 30, 2163–2175. Nesterov, Y., 2013. Gradient methods for minimizing composite functions. Mathematical Programming 140, 125–161. Parikh, N., Boyd, S., 2013. Proximal algorithms. Foundations and Trend in Optimization 1, 123–231. Patrinos, P., Guiggiani, A., Bemporad, A., 2015. A dual gradientprojection algorithm for model predictive control in fixed-point arithmetic. Automatica 55, 226–235. Rawlings, J.B., Mayne, D.Q., Diehl, M., 2017. Model predictive control: theory, computation, and design. 2nd edition ed., Nob Hill Publishing. Richter, S., 2012. Computational complexity certification of gradient methods for real-time model predictive control. Ph.D. thesis. ETH Zurich. Ryu, E.K., Boyd, S., 2016. A primer on monotone operator methods. Appl. comput. math. 15, 3–43. Saracco, P., Batic, M., Hoff, G., Pia, M.G., 2012. Uncertainty quantification (UQ) in generic MonteCarlo simulations, in: 2012 IEEE Nuclear Science Symposium and Medical Imaging Conference Record (NSS/MIC), IEEE. pp. 651–656. Simi´c, S., Bemporad, A., Inverso, O., Tribastone, M., 2022. Tight error analysis in fixed-point arithmetic. Formal Aspects of Computing 34, 1–32. Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S., 2020. OSQP: An operator splitting solver for quadratic programs. Mathematical Programming Computation 12, 637–672. Vakili, S., Langlois, J.M.P., Bois, G., 2013. Finite-precision error modeling using affine arithmetic, in: 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, IEEE. pp. 2591–2595. 8