{'fraction_non_alphanumeric': 0.0449301110944624, 'fraction_numerical': 0.013121561744118843, 'mean_word_length': 4.360027535566774, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Adaptive Partial Sensitivity Updating Scheme for Fast Nonlinear Model Predictive Control

In recent years, efficient optimization algorithms for Nonlinear Model Predictive Control (NMPC) have been proposed, that significantly reduce the on-line computational time. In particular, direct multiple shooting and Sequential Quadratic Programming (SQP) are used to efficiently solve Nonlinear Programming (NLP) problems arising from continuous-time NMPC applications. One of the computationally demanding steps for on-line optimization is the computation of sensitivities of the nonlinear dynamics at every sampling instant, especially for systems of large dimensions, strong stiffness, and when using long prediction horizons. In this paper, within the algorithmic framework of the Real-Time Iteration (RTI) scheme based on multiple shooting, an inexact sensitivity updating scheme is proposed, that performs a partial update of the Jacobian of the constraints in the NLP. Such update is triggered by using a Curvature-like Measure of Nonlinearity (CMoN), so that only sensitivities exhibiting highly nonlinear behaviour are updated, thus adapting to system operating conditions and possibly reducing the computational burden. An advanced tuning strategy for the updating scheme is provided to automatically determine the number of sensitivities being updated, with a guaranteed bounded error on the Quadratic Programming (QP) solution. Numerical and control performance of the scheme is evaluated by means of two simulation examples performed on a dedicated implementation. Local convergence analysis is also presented and a tunable convergence rate is proven, when applied to the SQP method.


I. INTRODUCTION
Nonlinear Model Predictive Control (NMPC) has been studied and applied intensively in the last decades. In NMPC, a nonlinear Optimal Control Problem (OCP) has to be solved on-line at every sampling instant. The OCP can be converted to a finite dimensional Nonlinear Programming (NLP) problem by direct methods, such as direct multiple shooting [1] and direct collocation [2]. The NLP problem can then be solved by a number of optimization algorithms, e.g., Interior Point Methods (IPM) [2] and Sequential Quadratic Programming (SQP) [3]. Fast NMPC algorithms based on direct methods have been proposed to speed up on-line optimization, see [4]- [6].
Efficient SQP algorithms based on direct multiple shooting for systems governed by Differential Algebraic Equations (DAE) have been thoroughly studied (see e.g. [7]). One of the computationally demanding steps of SQP methods when Yutao  applied to NMPC is the computation of sensitivities at each sampling instant, i.e. the Hessian of the Lagrangian and the Jacobian of the constraints. There are several methods for computing such sensitivities, e.g. finite difference [8], complex-step differentiation [9], and automatic differentiation [10].
Particularly, the Jacobian of constraints contains sensitivities of integration operators that parameterize continuous-time dynamics. Although efficient implementations of numerical integration with sensitivity generation are available [11], [12], sensitivity computation of this type still largely contributes to the overall on-line computational burden, especially for systems that are highly stiff or governed by implicit differential equations and DAEs.
In this paper, on one most promising SQP-based NMPC algorithms that is the Real-Time Iteration (RTI) [13], in which only one SQP iteration is performed at each sampling instant, is taken as the reference approach. The underlying idea is to initialize the new NLP by using information from the previous one, including states, controls, and multipliers making the closed-loop trajectory converging as system dynamics evolve, i.e. "on the fly" [14].
In the RTI framework with multiple shooting parameterization, a number of tailored approaches are available that employ suitable inexact sensitivities. In Multi-Level RTI (ML-RTI) [15], sensitivities are updated every m > 1 sampling instants. Hence, sensitivities are updated at a slower rate than other QP components. However, the choice of m is not intuitive and generally application dependent, thus requiring a long and complex tuning procedure. In ADJoint sensitivity RTI (ADJ-RTI) sensitivities computed off-line are used [16]- [19], and an adjoint sensitivity on-line computation is performed to identify the correct active-set and to ensure local convergence. Although computational cost is considerably reduced, thanks to the reduced number of sensitivity computations and condensing operations [18], this approach is effective only for systems exhibiting mild nonlinearities. Recently, partial sensitivity updating schemes called CMoN-RTI and DOPUS, that are tailored for multiple-shooting based NMPC, have been proposed [20]- [22]. In such schemes, the multi-stage feature of NLP problems arising in NMPC applications and the iterative nature of the solver are exploited. As a result, sensitivities are partially updated between two consecutive sampling instants. A so-called Curvature-like Measure of Nonlinearity (CMoN) or norm-criterion is used in a monitoring strategy to decide which and how many sensitivities should be updated. However, these monitoring strategies rely on heuristics and are strongly dependent on the application at hand.
In this paper, the partial sensitivity scheme CMoN-RTI of [21] is extended and improved. In particular, three main features are provided, namely: • a solution accuracy control strategy; • a practical tuning procedure; • convergence analysis.
From parametric optimization theory, the accuracy of the QP solution is related to parameters in the monitoring strategy. An advanced tuning strategy for CMoN-RTI is here developed that provides an automatic way to select which and how many sensitivities should be updated, while guaranteeing the QP solution a bounded Distance to Optimum (DtO) (i.e. the distance between the solution of the inexact sensitivity QP and that of the exact sensitivity QP). The tuning parameter is therefore the DtO tolerance, which has an important physical meaning.
The proposed scheme can significantly reduce the computational load when the system nonlinear dynamics are excited only on a small part of the prediction horizon (e.g., when regulating a system around its steady state or tracking a reference with look ahead). Moreover, since the additional computational time required by CMoN-RTI with respect to RTI is almost negligible, CMoN-RTI is a sensible alternative in all the scenarios where RTI is effective, as it usually yields an improvement in the average computational performance, hence saving computational power, and possibly an increase of the control frequency. In the worst case, CMoN-RTI degrades to RTI. A practical implementation of the scheme is given and its effectiveness is demonstrated by closed-loop simulations on two classical examples. An application of CMoN-RTI applied to the SQP framework with multiple iterations is also introduced. A tunable local convergence rate is proven.
The paper is organized as follows. In Section II, RTI and some inexact sensitivity schemes are briefly introduced to define the algorithmic framework. In Section III, the CMoN-RTI scheme is presented in detail. Section IV is devoted to the derivation of the advanced tuning strategy and to practical implementation aspects. In Section V, closed-loop simulation results using CMoN-RTI are shown. The CMoN-SQP is described in Section VI, and its convergence properties are discussed and demonstrated by a numerical example.

II. ALGORITHMIC FRAMEWORK
In this section, the standard RTI scheme [13] is introduced as the algorithmic framework of the paper. ML-RTI [15] and ADJ-RTI schemes [17] are here presented as two variants of RTI, with inexact sensitivity updating strategies.

A. Real-Time Iteration Scheme
In NMPC, a NLP problem can be formulated by applying direct multiple shooting [1] to an OCP over the prediction wherex 0 is the measurement of the current state. System states x k ∈ R nx are defined at the discrete time point t k for k = 0, . . . , N and the control inputs u k ∈ R nu are piece-wise constant. Here, (1d) is the inequality constraint where r(x k , u k ) : R nx × R nu → R nr . Equation (1c) refers to the continuity constraint where φ k (x k , u k ) is a numerical integration operator that solves the following initial value problem (IVP) 1 and return the solution at t k+1 .
The NLP problem (1) depends on the state and control initialization w and the state measurementx By embedding (1b) into (1c), the NLP problem can be written in a compact form as In RTI, problem (2) is solved by a tailored SQP method, where only one SQP iteration is performed at each sampling instant. At sampling instant i, the QP subproblem initialized at w i is defined as where ∆w = w − w i and ∇ is the gradient or Jacobian operator over w if no subscript is provided. The equality and inequality constraints are given by H i is the Hessian of the Lagrangian of (2), which is defined by L(w, λ, µ) := A(w)+λ ⊤ B(w)+µ ⊤ C(w), where λ, µ are Lagrangian multipliers associated with equality and inequality constraints, respectively. For most QP problems arising from NMPC, the Gauss-Newton Hessian approximation provides a sufficiently accurate Hessian with reduced computational burden [13]. Being independent of Lagrangian multipliers, the Gauss-Newton Hessian is adopted in this paper. Given the multi-stage nature of problem (1), matrices H i and ∇C(w i ) are block diagonal. In particular, the Jacobian matrix ∇B(w i ) has the following form: where ∇φ i k = ∇φ k (x i k , u i k ), k = 0, . . . , N − 1 and I a is an identity matrix of size a. The solution (∆w QP , λ QP , µ QP ) of (3) is used to update the solution of (2) by a single, full Newton step as Since only one QP problem is solved at each sampling instant, the RTI scheme is a special case of the linear, time-varying MPC strategy [14].

B. Inexact Sensitivity Schemes
When formulating problem (3), the Jacobian matrix ∇B(w) is computed at each sampling instant to obtain the current linearization of system dynamics [17]. Such computation involves sensitivity propagation of the numerical integration operator φ k in (1c) w.r.t. the initialization w k for each shooting interval, which can be computationally expensive for systems that are highly stiff or governed by implicit differential equations and DAEs. To avoid the repeated sensitivity computations, in ML-RTI, firstly proposed in [15], the computation of ∇B(w) is performed at a slower rate than other components, e.g. ∇A and B(w). Therefore, problem (3) with currently available but inexact sensitivities is solved at a faster rate [23], [24]. To account for the inexact Jacobian, a so-called optimality improvement step is employed by solving a slightly modified QP problem as whereb i = B(w i ) + ∇B∆w and ∇B is the Jacobian from the previous sample. The QP gradient is modified as where ∇Ã can be efficiently computed by applying adjoint sensitivity propagation schemes for ∇B ⊤ (w)λ [17], [18], which are much cheaper than the computation of the full Jacobian matrix ∇B(w). The multi-level framework with inexact sensitivities can be summarized in Algorithm 1 2 .
Algorithm 1 Multi-Level inexact sensitivity RTI scheme 1: Initialize (2) at (w 0 , λ 0 , µ 0 ). Choose a sensitivity update interval m ∈ N + . 2: for i = 0, 1, . . . do 3: 4: if i mod m = 0 then Solve (6) and obtain (∆w QP , λ QP , µ QP ) 9: Update the solution of the NLP problem by (5) 10: end for The ADJ-RTI scheme is a special variant of Algorithm 1 with m = ∞, in which the Jacobian matrix is computed only once off-line at the pre-defined initialization trajectory w 0 [17]. When applied in the SQP framework with multiple SQP iterations, such an inexact sensitivity scheme is proved to be convergent to the local minimum of the exact sensitivity NLP problem (2) [15]. The feasibility and stability of the adjoint scheme in SQP framework, without the optimality improvement, is analyzed in [19]. However, there remain some open issues when applying Algorithm 1 in RTI: • It is not trivial to choose an appropriate sensitivity update interval m, or a pre-defined trajectory w 0 , such that the inexact Jacobian ∇B is a good approximation of the exact one for every sampling instant i > 0. • The sensitivities over the prediction horizon, either updated or not, are treated as a whole. Therefore, the structure of the Jacobian matrix is not exploited. In the following Sections, a new sensitivity updating scheme is introduced that aims at overcoming the limitations described above.
III. CMON-RTI Several attempts to exploit the structure of the QP (3) are present in the literature. A Mixed-Level scheme has been proposed in [25], where only the first N c blocks in (4) are updated. However, choosing N c from heuristics may not be adequate for controlling highly nonlinear systems. Partial updating schemes where a fixed number of Jacobian blocks are updated have been independently proposed in [20] and [22] by using either CMoN or a so-called "norm criterion". In [21], an inexact scheme has been proposed, where a varying number of sensitivities are updated, namely, only the most "nonlinear" ones. In this paper, the CMoN-RTI scheme is extended by introducing adjoint CMoN on dual variables and analyzing the QP problem by using parametric optimization theory. In particular, an advanced tuning strategy is developed, that automatically provides the number of Jacobian blocks to be updated while granting that the DtO remains below a userdefined tolerance. For the sake of clarity, in the following the CMoN-RTI scheme proposed in [21] is summarized.

A. Curvature-like Measure of Nonlinearity
Studies on Measures of the Nonlinearity (MoNs) of nonlinear dynamic systems can be traced back to the 1980s. The three main classes of MoNs are: 1) the distance between a nonlinear system and its best linear approximation [26]; 2) the gap metric between two linear systems obtained by linearizing a nonlinear system around two different operating conditions [27]; 3) the curvature MoN (CMoN) at a point in the parameter space along a given direction. [28], [29]. Global and off-line metrics are developed in [26], [27]. CMoN is a local metric originally introduced to measure the nonlinearity in an estimation setting [28], [30], [31] and then extended to chemical processes control [29], [32]. It is defined as the ratio of the quadratic term over the linear term of the Taylor expansion of a nonlinear function z = g(s) along the ǫ direction in the input space: As the scaling effect of ǫ is cancelled out by using a square norm in the denominator of (7), this definition evaluates the instantaneous "curvature" of the manifold of z. However, a knowledge of up to second order derivatives of the function z, which are computationally expensive, is required. Also, higher order terms are not taken into account [33].
In [20], [21], a variant of CMoN has beeen proposed to measure the local nonlinearity of dynamic systems in the NMPC framework. Assuming that φ k in (1c) is twice differentiable in w k , the sensitivities of φ k w.r.t. the initialization at two consecutive sampling instants i and i − 1 satisfy where · denotes the Euclidean norm 3 and q i−1 is the distance between the two initializations. The tensor (8) is a vector of length n x with each element a (n x + n u ) by (n x + n u ) matrix. The computation of involves a vector-tensor-vector product and is defined in terms of n x vector-matrix-vector products [28]. The CMoN of φ k is defined by where higher order terms of φ k are included in the numerator.
Observe that, knowledge of only the first order derivative ∇φ i−1 k is required. According to (8), such CMoN measures the relative change of the directional sensitivities between two consecutive sampling instants. Observe that Similarly, an adjoint CMoN can be defined as follows to measure the relative change of the directional sensitivities over dual variables: 3 In the paper, all vector and matrix norms are Euclidean.
The term ∆λ i−1 ⊤ k+1 ∇φ i k can be computed by efficient adjoint sensitivity schemes. As will be shown in Sec. IV, (9) together with (10) play important roles in controlling the accuracy of QP solutions.
At each sampling instant, the nonlinearity of a dynamic system over the entire prediction horizon can be estimated by applying (9) and (10) to each shooting interval.

B. Updating Logic
Due to the multiple shooting discretization, each block ∇φ i k in (4), and the corresponding CMoN κ i k , uniquely depend on the initialization w k . Hence, evaluation of CMoN, integration, and sensitivity generation can be performed independently at each shooting interval . Thus, the set of sensitivity blocks {∇φ i k } can be divided into two parts: 1) an updating subset where the sensitivity blocks are updated; and 2) the remaining subset where the sensitivity blocks are kept unchanged. If the first subset is much smaller than the second one, a significant reduction of computational cost for sensitivity evaluations can be achieved. To this end, CMoN can be used to determine such an updating subset. Intuitively, when κ i k is sufficiently small, the sensitivity ∇φ i k is close enough to ∇φ i−1 k , hence sensitivity update is not necessary for the current sampling instant. The block k of the Jacobian matrix ∇B is updated according to the following strategy. Set the values of thresholds η i pri and η i dual , where the subscript pri denotes the primal variable and dual the dual variable. Then, The proposed strategy is effective in both of the following cases: 1) Regulation: Given a sufficiently long prediction horizon, the system nonlinearity is typically excited in a small part of the predicted trajectory, that is, far from the steady state. As the system is approaching its steady state, less and less sensitivities are expected to be updated. 2) Reference tracking: Assuming that future reference is known in advance, a widely used choice is to progressively update the reference starting from the end of the prediction horizon [14]. This approach has the beneficial impact that the initial part of the predicted trajectory is not affected by the reference change. As a result, given a sufficiently high sampling frequency, sensitivity update is necessary only in the final part of the prediction horizon, whereas information from the past can be effectively used elsewhere.

IV. AN ADVANCED TUNING STRATEGY
In (11), thresholds η pri and η dual regulate the trade-off between the accuracy of the Jacobian approximation and the computational cost, by determining the updating subset with the largest CMoN values. An intuitive way of choosing the thresholds is to set a constant value, i.e. η i pri = η 0 pri and η i dual = η 0 dual , for all sampling instants. When η 0 pri = 0 and η 0 dual = 0, the proposed scheme becomes the standard RTI scheme with N f = N , i.e. all sensitivities are updated at every sampling instant. When η 0 pri ≥ max(κ i k ) and η 0 dual ≥ max(κ i k ) for all i, N f = 0 and no sensitivity is updated on-line, hence CMoN-RTI coincides with ADJ-RTI [15]. Thresholds η 0 pri and η 0 dual can take any value in the sets [0, max(κ i k )] and [0, max(κ i k )], respectively, to achieve a flexible tuning. A tuning strategy suitable for real-time implementation can be used: η i pri and η i dual can be chosen to update, at each instant, a fixed number of sensitivities [20], [22]. However, a pre-defined limited number of sensitivity updates may not be suitable for controlling highly nonlinear systems.
A satisfactory trade-off between the accuracy of the sensitivity approximation and computational cost can be achieved by means of an advanced, time-varying tuning of the thresholds η pri and η dual . The key observation is that using the inexact Jacobian in (6) affects the accuracy of both primal and dual solutions. A relation that reflects inaccuracy of the sensitivities into inaccuracy of the solution of the QP problem can therefore be used to choose, at each sampling instant, the values of η pri and η dual , that guarantee a tunable, bounded error on the QP solution. By adopting this strategy, CMoN-RTI can adjust the number of updated Jacobian blocks according to system operating conditions to achieve a numerical and control performance as close as possible to the standard RTI scheme, with improved computational performance.
First some facts from parametric programming theory are reviewed, then the advanced tuning strategy is detailed and some practical implementation aspects are finally considered.

A. Parametric Nonlinear Programming: stability of the solution
Two definitions concerning parametric QP are first introduced. The Jacobian approximation error is taken as a perturbation parameter. Three Lemmas describing the stability of the QP solution w.r.t. to such parameter are then given. Definition 1. Define a parametric QP(p) with parameter vector p ∈ R np in the equality constraint as where ∇L is the gradient of the Lagrangian of (2), b(p) = B(w) + (∇B + P )∆w, P := ∇B − ∇B is the Jacobian approximation error, and ∇B is the inexact Jacobian with partially updated blocks. The perturbation vector p = vec(P ) ∈ R np is the vectorization of P after eliminating zero elements.
According to Definition 1, the exact Jacobian QP problem (6) is referred to as QP(0). Due to multiple shooting discretization, P has the following banded block structure where O a is a zero matrix of dimension a and P k ∈ R nx×(nx+nu) is the k−th block of the Jacobian approximation, and in general is a dense matrix.
Observe that QP (13) has a modified objective gradient with respect to (6). However, it can be easily proved that these two formulations are equivalent [34]. The additional computational cost can be neglected since both formulation (13) and (6) contain adjoint sensitivities in their objective. We adopt (13) as the standard form hereafter.
The following Lemma shows that the distance between the primal solutions of QP(0) and QP(p) is bounded, and the bound is of the same order of the Jacobian approximation error. Lemma 1. [35]. Let ∆w(0) and ∆w(p) minimize QP(0) and QP(p) over corresponding feasible sets, respectively. Then there exists constants c and ǫ * > 0 such that ∆w(p) − ∆w(0) ≤ cǫ whenever ǫ ≤ ǫ * and ǫ = ∇G − ∇G = P .
The following Lemma shows that the solution ∆y(p) is a unique minimizer of (13). Moreover, the active set is locally stable. Lemma 2. [36] Under the assumption on differentiability, second-order sufficient conditions, constraints linear independence and the strict complementary slackness condition, there exists a unique solution ∆y(p), which is continuously differentiable w.r.t. p for p in a neighborhood of 0. Moreover, the set of active inequality constraints is unchanged, strict complementary slackness holds, and the active constraint gradients are linearly independent at ∆w(p).
Finally, the following Lemma provides a linearly approximated relationship between the exact and inexact solutions. where M, N are given in (12), and b k and c k are the k−th row of b(p) and c, respectively. Lemma 3. Modern studies based on perturbation theory show that the solution manifold ∆y(p) is nonsmooth but continuous. Therefore, ∆y(p) can be close enough to ∆y(0) even in the presence of active-set changes. The reader is referred to [37], [38] and references therein for more details.

B. First Order Error Analysis
In the neighborhood of p = 0, (3) can be rewritten as As shown in Appendix A, it holds that and, by pre-multiplying (14) by M (p), it follows that M (p)∆y(0) = M (p)∆y(p) + N (p)p.
Therefore, the DtO at the sampling instant i satisfies Note that, given a finite dimensional and non-singular real matrix M (p i ), its Euclidean norm is bounded. Hence, the DtO e i is bounded only if P i ⊤ ∆λ(p i ) and P i ∆w(p i ) are bounded. The two bounds are referred as the dual bound and primal bound respectively, and are discussed in the following. 1) Primal Bound: By using the primal threshold η pri in the updating logic (11), one obtains where pri is a vector of directional sensitivities given by Hence, a bound in the direction of the primal variable is as follows 2) Dual Bound: Similarly, for adjoint CMoN it holds that Hence, a bound in the direction of the dual variable is obtained as follows

C. Thresholds Estimation
Given a DtO toleranceē i at the sampling instant i, let where 0 < c 1 < 1 is a tuning parameter that trades off impact of the primal and dual bounds on the DtO. By substituting (20) into (15), one obtains e i 2 ≤ ē i 2 . Therefore, the primal and dual thresholds satisfy the following inequalities: Theorem 1. U 1 , U 2 : R × R → R are piecewise discontinuous functions of (η i pri , η i dual ) and their ranges are finite sets. The proof of Theorem 1 is given in Appendix C. Given Theorem 1, a formal solution to find the maximal (η i pri , η i dual ) is then to solve the following problem The solution of problem (22) provides the maximal values of the thresholds, corresponding to the minimum number of sensitivity updates while guaranteeing a bounded DtO. Note that for a givenē i ≥ 0, there always exists at least one feasible solution to (22), i.e. (η i pri , η i dual ) = 0, that makes CMoN-RTI coincide with the standard RTI scheme.

D. Practical Implementation
Problem (22) can be solved via enumeration, which requires to repeatedly solve problem (13). However, this is computationally prohibitive and undermining the advantage of CMoN-RTI. A practical approach to avoid solving problem (22) is setting the two thresholds at their upper bounds in (21), by using approximated information from previous sampling instants.
Firstly, the unknown ρ i in (16) is replaced by ρ 0 . The rationale for such choice is given by the fact that ρ i is the reciprocal of the smallest singular value of M (p i ). According to (12), M (p i ) is a very sparse matrix and its smallest singular value is close to 0 and does not vary much between sampling instants. Hence, ρ 0 cab be computed offline and used for all on-line computations.
Secondly, as shown in (18) and (19), the values of (α i , β i ) cannot be computed in a real time implementation, since the Jacobian approximation error P i cannot be computed from approximate sensitivities, and the solution (∆w(p i ), ∆λ(p i )) is not known in advance. When the NMPC controller is converging on the fly, it holds that ∆w(p i ) ≤ q i−1 . In such case, α i and β i are usually less than one. Since larger values of α i and β i give more conservative results (as they lead to the computation of a larger number of sensitivities), a sensible choice is setting (α i , β i ) = (1, 1). This aspect is also discussed in Section V and VI with reference to practical implementation of the algorithm.
Thirdly, c 1 is the parameter that allows to balance the impact of the primal and dual thresholds on the DtO. In this paper, the choice c 1 = 0.1 is made since the magnitude of multipliers is typically bigger than the primal solution. As shown in Section V, this choice allows to achieve a satisfactory performance.
Fourthly, in (15), an upper bound for DtO is obtained by means of norm inequality. Such inequality may lead to conservative upper bounds of thresholds (21), hence updating more sensitivities than necessary. To account for such issue, a scaling parameter γ i is introduced such that The value of γ i cannot be computed on-line by using (23) for real time applications. However, an estimate of it can be obtained by relying on the following theorem.
Theorem 2. For a real, linear system z = Xt where z ∈ R m , t ∈ R m , X ∈ R m×m , it holds that z = X · t if and only if Σ = σ 2 I, where X = U ΣV ⊤ is the Singular Value Decomposition (SVD) of X.
Theorem 2 can be proved by applying the definitions of SVD and spectral norm of matrices. According to Theorem 2, if the singular values of M −1 (p i ) are all equal, (23) holds for γ i = 1. For general matrices whose singular values are not identical, (23) holds for γ i > 1. Hence, γ i is estimated by where std(Σ) is the standard deviation operation and σ(M −1 (p i )) is the set of singular values of M −1 (p i ). To make on-line computation feasible, γ 0 can be used, as it can be computed off-line. Effectiveness of this choice is discussed in Section V. Finally, the approximated thresholds estimates are given by A summary of the practical implementation of CMoN-RTI is given in Algorithm 2.

V. NMPC SIMULATION CASE STUDY
In this section, Algorithm 2 is applied to two examples, namely, the control of an inverted pendulum and of a chain of masses. Numerical integration and sensitivity generation are performed by a 4 th order explicit Runge-Kutta integrator with 4 steps per shooting interval, provided by the CasADi toolbox [39] using automatic differentiation. The QP problem is solved by using HPIPM, a structure-exploiting interior point solver based on hardware tailored linear algebra libraries [40]. Algorithmic parameters are chosen as described in Section IV.D for all examples. The computing environment is Ubuntu 16.04 on a PC with Intel core i7-4790 running at 3.60GHz, and the implementation is coded in plain C with -O2 compilation optimization flag.

A. Inverted Pendulum
An inverted pendulum is mounted on top of a cart and can roll up to 360 degrees. The dynamic model is given bÿ p = −m1l sin(θ)θ 2 + m1g cos(θ) sin(θ) + F m2 + m1 − m1(cos(θ)) 2 , θ = 1 l(m2 + m1 − m1(cos(θ)) 2 ) (F cos(θ) − m1l cos(θ) sin(θ)θ 2 + (m2 + m1)g sin(θ)), where p, θ are the cart position and swinging angle, respectively, and F is the control force acting on the cart. The model and values of parameters m 1 , m 2 , l, g are taken from [41]. For this example, a time-varying reference is given to the inverted pendulum to track different horizontal displacements and swing angles. A perfect initialization is chosen by optimally solving the OCP for t = 0 off-line. A short (N = 40) and a long (N = 120) prediction horizon are applied with a control interval T s = 0.05s. The tolerance on DtO in CMoN-RTI follows the rule given bȳ where (ǫ abs , ǫ rel ) are the absolute and relative tolerances, n is the number of optimization variables, and y = (x, λ, µ) is the optimal triple. Such choice ensures that the DtO tolerance scales with the size of the problem and the scale of the variable values [42]. For this problem, we set ǫ abs = 10 −1 , ǫ rel = 10 −1 . In Figure 1, the closed-loop state and control trajectories generated by the standard RTI and CMoN-RTI with the two prediction horizons are shown. The control performance of CMoN-RTI is indistinguishable from that of the standard RTI scheme, which demonstrates that CMoN-RTI is able to maintain the closed-loop performance as the standard RTI while using much less sensitivity computations. In Figure 2, the percentage of exactly computed sensitivities per sampling instant is given. The CMoN-RTI scheme can adapt to operating conditions by evaluating more sensitivities when the reference is about to change, as the peaks occur at around t = 3, 8, 13, 18s. A significant reduction of the percentage of updated sensitivities is observed when N = 120, making CMoN-RTI adequate to deal with the case of long prediction horizons. As explained in Section III-B , only the last part of the reference is triggering sensitivity updates. Hence, the longer the prediction horizon, the lower the percentage of sensitivities to be updated. Figure 3 shows the DtO at each sampling instant together with the user-defined tolerance. An additional QP with exact Jacobian matrix is solved at each sampling instant to compute the DtO. In both cases, the DtO is lower than the tolerance.
To examine the effectiveness of using ρ 0 and γ 0 in (25), the relative difference is defined as For the inverted pendulum example, it is observed that the maximum value of r is 24%, i.e. a sufficiently small difference  (Fig. 4a) and the control trajectories (Fig. 4b) in one of the simulations using the standard RTI scheme. One end of the chain is fixed on a wall, while the other end is free and under control. The control interval is Tc = 0.2s. The control inputs are constrained by u(t) ∞ ≤ 1.  Fig. 5. The time tst needed to stabilize the chain of masses using RTI, ML-and CMoN-RTI in the total 50 simulations. The chain of masses is considered to be stabilized at time tst that is computed by (29). For all schemes, the stabilizing time is set to be tst = 50s if the chain is not stabilized within 50s.
that confirms the effectiveness of the approximation strategies discussed in Section IV.

B. Chain of Masses with Nonlinear Springs
A chain of masses is a system with n masses connected by springs on a chain [16]. The dynamic model is given bẏ where p i (t) ∈ R 3 and v i (t) ∈ R 3 are the positions and velocities of the i−th mass, respectively, and is the spring force from mass i to i + 1 and F N L is its nonlinear component. The velocities of the free massṗ n (t) are controlled by u(t). As demonstrated in [16], [18], ADJ-RTI is able to stabilize the chain of masses if F N L = 0, i.e. the chain is connected by linear springs. In this paper, nonlinear springs [43] are considered with A total of 50 simulations are performed while using the standard, ML-, ADJ-and CMoN-RTI, with randomly assigned initial positions and velocities of the masses, see e.g. Fig 4, for the positions and control trajectories generated in one of the simulations. For ML-RTI, the entire constraint Jacobian matrix is updated every m = 2 sampling instants; For ADJ-RTI, the Jacobian matrix is computed off-line at the steady state trajectory; For CMoN-RTI, the DtO tolerance is chosen as in (28). To ensure that an accurate representation of the system is always used in the controller, at least 10% sensitivities are updated at each sampling instant. These sensitivities are those having the largest values of CMoN, hence exhibiting the most significant nonlinearities [20]. Control performance, numerical robustness, and efficiency of CMoN-RTI are evaluated and compared with standard RTI, ML-RTI, and ADJ-RTI. Firstly we collect statistics of the stabilizing time t st , defined as from the 50 simulations. In Fig. 5, the statistics of the standard, ML-, ADJ-and CMoN-RTI with N = 40, 80, 160 are shown. Note that, if the chain is not stabilized within 50s, we set t st = 50s, which is a conservative choice since the stabilization process may take far more than 50s. For all simulations, RTI is able to stabilize the chain within 50s.   The mean and interquartile range (IQR) of t st of CMoN-RTI is very close to those of the standard RTI. This means that CMoN-RTI has a similar control performance to the standard RTI in most of the situations. On the other hand, ML-RTI has a similar stabilizing time to RTI when N is short, whereas t st grows significantly as N becomes larger. ADJ-RTI, initialized at the steady state trajectory, is not able to provide acceptable control performance, especially when N is large. The control performance is also evaluated by assessing the optimality of each controller. In Fig. 6, the average Karush-Kuhn-Tucker (KKT) value, i.e. the norm of the gradient of the Lagrangian of the NLP (2), at each sampling instant is presented. It can be observed that • ML-RTI KKT values exhibit strongly oscillatory behavior since the Jacobian update is performed every m = 2 sampling instants only. • As the system converges "on the fly", the KKT of CMoN-RTI decreases smoothly as that of the standard RTI. As for numerical robustness, the number of simulations where each controller fails to stabilize the system within 50s is reported in Table II . Given that the initial condition of each simulation is randomly assigned, the numerical robustness or the sensitivity w.r.t. initialization of each controller can be assessed. CMoN-RTI is able to stabilize the chain within 50s in all situations, although the maximal stabilizing time is larger than that of the standard RTI (see Fig. 5). ML-RTI has a few failed cases if N = 40 and this number increases as the prediction horizon grows. Not surprisingly, ADJ-RTI exhibits the poorest robustness properties as it heavily depends on the quality of the off-line Jacobian matrix.
To evaluate computational efficiency, the average percentage of exactly updated sensitivities using CMoN-RTI at every sampling instant is reported in Figure 7. In the first 20s, the KKT values of CMoN-RTI and RTI are almost identical, however, the number of updated sensitivities is at most 80% and it reduces to 60% when N becomes larger. After t = 30s, when the system is close to its steady state, updating only 10% of the blocks allows to still maintain small KKT values. Table I shows the average and maximal computational time of CMoN-RTI and the RTI scheme per sampling instant. For this example, the speedup factor is computed by using the maximal computational time. With different prediction horizons, the computational time for evaluating sensitivities varies from about 60% to 36% of the full RTI step, and the one for the QP varies from about 36% to 65%. As a result, the speed up factor is strongly related to the distribution of computational time among the critical steps of the full RTI. In the specific example, the speedup factor is always positive with a maximum value of 14.2% when N = 80. Also, observe that the computational performance obtained in the examples is related to the use of the simple explicit Runge-Kutta integrator. For systems that require the use of more complex integrators, whose sensitivities are more computational expensive, CMoN-RTI is expected to achieve a greater speedup factor.
VI. CONVERGENCE ANALYSIS Algorithm 2 is a partial sensitivity updating scheme in the framework of RTI between two consecutive sampling instants. It can also be straightforwardly extended to the SQP framework, where a sequence of QP problems is solved until convergence is achieved. The resulting algorithm, denoted hereafter as CMoN-SQP, partially updates sensitivities between two consecutive SQP iterations. In the SQP scenario, the Two-side-Rank-One (TR1) updating SQP algorithm has been proposed in [44] for equality constrained problems. Similar to the famous Symmetric-Rank-One (SR1) updating scheme [8], the TR1 scheme requires Hessian and Jacobian updates to satisfy both direct and adjoint secant conditions. This method is extended to linearly inequality constrained problems in [34] and its local convergence is proved.
Differently from the TR1 scheme, which adopts a rank one Jacobian matrix update, CMoN-SQP achieves a block update by exploiting the structure of the problem. In addition, the primal and dual bounds are satisfied, instead of enforcing secant conditions. In the following, local convergence of CMoN-SQP is proved and it is shown that the convergence rate is tunable via the choice of the DtO tolerance.

A. Local Convergence of CMoN-SQP
Consider the parametric QP problem (13). Solving problem (13) in a SQP algorithm is equivalent to solving the following nonlinear system: where λ denotes the multiplier for both equality and active inequality constraints, C a contains the active constraints, and R is a matrix with orthonormal column vectors, such that ∇C a R = 0 [34]. The Jacobian matrix of the nonlinear system is where H i is an approximation of the exact Hessian, i.e. the Gauss-Newton approximation which is independent of the multiplier λ. Let J i be an approximation of the exact Jacobian ∇F (y i ) with The following theorem indicates that the proposed scheme is convergent in the neighborhood of p = 0.
Theorem 3. Let F : V → R ny , V ⊂ R ny be continuously differentiable. Consider the two sequences Assume that 1) the Jacobian matrix is invertible, uniformly bounded, and has uniformly bounded inverses, 2) there exists a κ 0 < 1 such that ∆y i+1 * ≤ κ 0 ∆y i * for all i > m 1 , m 1 ∈ N. Hence, starting from y 0 ∈ V, the sequence {y * } converges to a local optimizer y + * , 3) J(y i p ) is generated by Algorithm 2, Then,

1) there always exists a set of scalars
such that the distance between the sequences {y p } and {y * } is sufficiently small at each iteration, 2) there always exists a set of scalars {i ∈ N + |ē i ≥ 0} and a κ 2 satisfying κ 0 ≤ κ 2 < 1, such that ∆y i+1 p ≤ κ 2 ∆y i p for all i > m 2 , m 2 ∈ N, and the sequence {y p } converges to y + p = y + * starting from y 0 . Proof: Let the locally exact solution initialized at y i p be Assume that at iteration i, the DtO is satisfied as Let d i y = y i p − y i * be the distance between the two sequences at the current iteration. Observe that Assume that d i y is sufficiently small and O( d i y 2 ) can be neglected, then by combining (30) and (31), it follows that As a result, where g i = ∇F −1 (y i * )(∆y i ⊤ 0 ∇ 2 F (y i * )) + I . The distance between the two solutions at the current iteration is and the distance between the two sequences at the next iteration is Since Algorithm 2 always starts from d 0 y = 0, d i y , ∀i > 0 is a linear combination of (ē 0 ,ē 1 , . . . ,ē i ). Therefore, it is always possible to choose a set of scalars {i ∈ N + |ē i ≥ 0}, such that d i y ≈ 0. Equivalently, the sequence {y p } can be sufficiently close to {y * } at every iteration.
Consider now the convergence properties of {y p }. By assumption 2, it follows that can be arbitrarily small, there exists a κ 2 satisfying κ 0 ≤ κ 2 < 1 such that Therefore, the sequence {y p } is convergent and its convergence rate is at most identical to that of {y * }. As proved in [15], [34], when {y p } does converge, it converges to the exact limit y + * of the sequence {y * }. Theorem 3 shows that the Jacobian approximation error can be controlled by using user-defined DtO tolerances, hence the convergence can be satisfied by using appropriate tuning configurations. The convergence rate is also shown to be tunable, which increases the flexibility of the proposed algorithm. Ifē i = 0, ∀i ≥ 0, CMoN-SQP becomes the standard SQP algorithm with the same convergence rate.

B. Numerical Examples
As an example, the CMoN-SQP scheme is applied to the inverted pendulum (26). The control objective is to invert the pendulum from bottom to top. CMoN-SQP is used to solve the OCP in open-loop at time t = 0 with N = 40. Since only local convergence is of interest, the initialization of the OCP is in a neighborhood of the optimal solution and a full Newton-step is adopted at each iteration. Figure 8 shows the convergence behavior of two different DtO choices of CMoN-SQP. The left y-axis reports the KKT value, that indicates the optimality of the solution. The right y-axis reports the percentage of sensitivities being updated at each iteration. To show how the choice of DtO tolerance affect the convergence rate, the following DtO tolerances are used: s 1 : (ǫ abs = 10 −2 , ǫ rel = 10 −2 ) s 2 : (ǫ abs = 10 −1 , ǫ rel = 10 −1 ) A more aggressive setting (s 2 ) leads to less sensitivity evaluations but slower convergence rate. Figure 9 shows the convergence behavior of three choices of c 1 for (25). The convergence rate is not sensitive to the values of c 1 . In practice, to solve a structured NLP problem by using CMoN-SQP, one would achieve a satisfactory trade-off between the cost of sensitivities and the number of iterations by properly tuning the DtO tolerance.

VII. CONCLUSION
In this paper, the partial sensitivity updating scheme CMoN-RTI of [21] is extended by proposing an advanced tuning strategy with solution accuracy control and convergence analysis. In CMoN-RTI, sensitivities are updated based on CMoN of the dynamic system over the prediction horizon. The CMoN works as a metric to evaluate the quality of sensitivity approximation, and an updating logic is designed with the use of primal and dual thresholds. By using parametric optimization theory, an advanced strategy for tuning such thresholds is proposed. Such strategy automatically determines the number of updated sensitivities and guarantees the DtO of QP solutions under a user-defined tolerance.
Closed-loop simulations show that the CMoN-RTI scheme exhibits good control performance when applied to highly nonlinear control problems. The application to an inverted pendulum shows that CMoN-RTI can adapt to reference changes while satisfying the DtO tolerance. The results from a chain of masses with nonlinear springs demonstrate the superior control performance, numerical robustness, and efficiency of CMoN-RTI.
The proposed scheme has also been extended to full SQP algorithms, denoted as CMoN-SQP and its local convergence is proved. Comparing to existing inexact sensitivity SQP methods, CMoN-SQP has two unique properties, namely, tunable convergence rate and structure exploiting updating logic.
Future studies may focus on possible extensions of CMoN-RTI. While ADJ-RTI and Mixed-Level schemes can benefit from condensing steps with significantly reduced computational efforts [18], [25], CMoN-RTI usually requires to perform a full condensing step at every sampling instant. Further improvements can be achieved by adopting partial condensing methods [45]. where P j,: is the jth row of P . For elements in N , it holds where O 1 ∈ R (nx+nu)×(nx +nu) , O 2 ∈ R (nx+nu)×nx , Λ j = I nx+nu ⊗ ∆λ j and W k = In x ⊗ ∆w ⊤ k with ∆w ⊤ k ∈ R 1×(nx +nu) .
APPENDIX B In this Appendix, the expression for V i−1 pri in (17) is derived. From the updating logic (11), it can be easily obtained that For the full Jacobian matrix, it holds A similar derivation can be conducted for V i−1 dual . The details are hence omitted.

APPENDIX C
A proof of Theorem 1 is given in this Appendix.
It follows that sensitivities with larger CMoN values always get updated first. As a consequence, there are finite number of combinations (actually N + 1) of possible updated sensitivities, resulting in a finite number of possible P i matrices. Therefore, ρ i , α i , β i which are defined in (16), (18) and (19), are functions of P i and have at most N + 1 possible values. Hence, the range of U 1 , U 2 are finite. In addition, for some k = j, ∀(η i pri , η i dual ) satisfying there exist (ǫ pri , ǫ dual ) = 0 ∈ R such that κ i k < η i pri + ǫ pri < κ i j , κ i k < η i dual + ǫ dual <κ i j .
Hence, P i remains constant under the perturbation of (ǫ pri , ǫ dual ). The values of U 1 , U 2 are constant for any (η i pri , η i dual ) satisfying (32). Discontinuity exists at (η i pri , η i dual ) = (κ i j ,κ i j ), when the matrix P i has different number of nonzero blocks resulting in different values of U 1 , U 2 .