Complex Eng SystComplex Engineering SystemsOAE Publishing Inc.10.20517/ces.2021.09Research ArticleA qLPV Nonlinear Model Predictive Control with Moving Horizon EstimationMoratoMarcelo Menezes^{1}^{2}http://orcid.org/0000000271370522BernardiEmanuel^{3}StojanovicVladimir^{4}^{1}Departamento de Automação e Sistemas, Universidade Federal de Santa Catarina, Florianópolis 88040900, Brazil.^{2}Univ. GrenobleAlpes, CNRS, Grenoble INP (Institute of Engineering, Univ. GrenobleAlpes), GIPSALab, Grenoble 38000, France.^{3}Applied Control & Embedded Systems  Research Group (AC&ESRG), Universidad Tecnológica Nacional, San Francisco 2400, Argentina.^{4}The Faculty of Mechanical and Civil Engineering in Kraljevo, Department of Automatic Control, Robotics and Fluid Technique, University of Kragujevac, Dositejeva 19, Kraljevo 36000, Serbia.Correspondence Address: Marcelo Menezes Morato, DAS/CTC/UFSC, Trindade, Caixa Postal 476, FlorianópolisSC 88040900, Brazil. Email: marcelomnzm@gmail.com.
Received: 27 Aug 2021  Revised: 22 Sep 2021  Accepted: 28 Sep 2021  Published: 15 Oct 2021
This paper presents a Model Predictive Control (MPC) algorithm for Nonlinear systems represented through quasiLinear Parameter Varying (qLPV) embeddings. Inputtostate stability is ensured through parameterdependent terminal ingredients, computed offline via Linear Matrix Inequalities. The online operation comprises three consecutive Quadratic Programs (QPs) and, thus, is computationally efficient and able to run in realtime for a variety of applications. These QPs stand for the control optimization (MPC) and a MovingHorizon Estimation (MHE) scheme that predicts the behaviour of the scheduling parameters along the future horizon. The method is practical and simple to implement. Its effectiveness is assessed through a benchmark example (a CSTR system).
Model Predictive Control (MPC) is a very powerful control method, with widespread industrial application. The core idea of MPC ^{[1]} is simple enough: a process model is used to predict the future output response of the process; then, at each instant, the control law is found through the solution of an online optimization problem, which is written in terms of the model, the process constraints and the performance goals. For the case of processes represented by Linear TimeInvariant (LTI) models, MPC is translated as a constrained Quadratic Programming Problem (QP), which can be evaluated in realtime by the majority of standard solvers.
Extra attention should be payed to the fact that the theoretical establishment MPC was especially consolidated after the proposition of “terminal ingredients”, which served to demonstrate robust stability and recursive feasibility properties ^{[2]}. These properties are enabled when some conditions with respect to a terminal stage cost V(·) and to a terminal constraint X_{f} are verified. Essentially, the terminal set must be robust positively invariant for the controlled system, the stage cost must be class lower bounded and the terminal cost V(·) should be class upper bounded and Lyapunovdecreasing (it must decay along the horizon).
For many years, MPC was mostly seen in the process industry, regulating usually slower applications (with longer sampling periods). This was mainly due to the fact that the inherent optimization procedures were excessively costly (numericalwise) and became impractical for realtime systems.
Nonlinear MPC (NMPC) algorithms yield complex optimization procedure, with exponential growth of the numerical burden. Nevertheless, the majority of system is indeed nonlinear and, thus, literature has devoted special attention to feasible NMPC design since the 00’s ^{[3]}. Originally, NMPC algorithms were hardly able to run in realtime ^{[4]}, but recent research effort has focused to a great extent on ways to simplify or approximate, usually through GaussNewton, Lagrangian or multipleshooting discretization approaches ^{[5]}, the online Nonlinear Programming Problem (NP) in order to make it viable for fast, timecritical processes. Some of these faster NMPC algorithms run within the range of a few milliseconds, resorting to solverbased solutions (as in ACADO^{[6]} or GRAMPC^{[7]} algorithms) or GPUbased schemes^{[8,9]}.
Parallel to these approximated methods, another research route is now expanding to address the complexity drawback of “fullblown” NMPC strategies: using quasi/Linear Parameter Varying (qLPV/LPV) model structures to embed the nonlinear dynamics, as in ^{[10]}, and thus facilitate the online optimization. Since LPV models retain linearity properties through the input/output channels, the optimization can be reduced to the complexity of a QP. A recent survey^{[11]} details the vast possibilities of issuing NMPC through LPV structures. The basic requirement of these methods is that the nonlinearities must respect the Linear Differential Inclusion (LDI) property^{[12,13]}, in such a way that they can be embedded into a qLPV realisation, appropriately “hidden” in scheduling parameters ρ.
Instead of using a movingwindow linearization strategy to yield fast NMPCs with timevarying models^{[14]}, or of using approximated solutions of the NP iterations ^{[6]}, this paper follows the lines of the qLPV embedding framework, which allows for an exact description of the nonlinear system and, thereby, no timeconsuming linearization or Jacobian computation needs to take place. As previously evidenced^{[11]}, these qLPV methods are able to use the scheduling proxy ρ(k) = f(x(k), u(k)) to compute the process predictions rapidly. In fact, these methods have recently been shown ^{[15]} to outrank (or perform equivalently as) fast NMPC solvers, such as ACADO. Some of these recent development are further detailed:
Some works^{[16,17]} opt to consider a frozen/constant guess for the scheduling parameters along the future horizon and ensure, through the use of terminal ingredients, that the trajectories are sufficiently regulated, despite the uncertainty along the horizon;
Morato et al. ^{[18]} propose a method to determine an educated estimation for scheduling variables using a recursive LeastSquares procedure. A similar procedure is applied in^{[19]}. The main drawback is that the results could be suboptimal, meaning that local minima found through their QPs/ Sequential QPs (SQPs), which may not ensure sufficient performances.
The most prominent results are those reported in the recent works by Cisneros & Werner ^{[2022]}. The original idea^{[20]} is to iteratively use the prediction for the future state trajectories (output of the QP) to compute a guess of the scheduling parameters, using the nonlinear proxy ρ(k + j) = f(x(k + j)). The method was extended ^{[21]} to referencetracking and shown to yield a Secondorder Cone Program (2^{nd}OCP) formulation for the resulting NMPC, which is easier to solve than an NP. The formulation was further smoothed in the most novel reference ^{[22]}, wherein the procedure is split into an offline preparation part, using Linear Matrix Inequalities (LMIs) to compute a robust positively invariant terminal set, and an online residing solely in reiterating SQPs.
1.1. Contributions and rganization
oAs detailed in the prequel, the topic of NMPC through qLPV embedding has been studied by a handful of papers and deserves further attention. It seems that the development of these strategies can surely be established as a competitive category for nonlinear MPC design, regarding timecritical applications.
Pursuing this matter and motivated by the previous discussion, this paper proposes an alternative formulation to the recent algorithm by Cisneros and Werner^{[22]}. In their work, the nonlinear proxy has to be evaluated online w.r.t. to the future state evolution prediction originated through the QP. The alternative procedure proposed herein relies on approximating the nonlinear proxy by a timevarying autoregressive function, whose parameters are found through another QP, based on a MovingHorizon Estimation (MHE) method. This alternative is able to slightly boost the numerical performances of the whole algorithm, which only needs to evaluate three QPs to find the control law.
Accordingly, the contributions presented are the following:
An alternative formulation for NMPC is proposed: using qLPV embeddings, the MPC operates together with an MHE layer, which estimates the future behaviour of the scheduling parameters.
The convergence of the algorithm is demonstrated.
A benchmark example is used to demonstrate the effectiveness of the proposed scheme, in terms of performance and numerical burden.
Regarding organization, this paper is structured as follows. Idn the next Section, the preliminaries and formalities are presented, especially regarding how nonlinear processes can be embedded into a qLPV representation through LDI. Moreover, the problem setup regarding MPC applied to such qLPV model is presented. Furthermore, the proposed MHEMPC formulation and the discussion about stability and an offline LMIsolvable remedy for the computation of the terminal ingredients is addressed. Lastly, Section simulation results and general conclusions are drawn.
1.2. Basic efinitions
Definition 1. Nonlinear Programming Problem
Consider an arbitrary realvalued nonlinear function f_{c}(x_{c}). A Nonlinear Programming Problem determines the vector x_{c} that minimizes f_{c}(x_{c}) subject to g_{i}(x_{c}) ≤ 0, h_{j}(x_{c}) = 0 and , where g_{i} and h_{j} are also nonlinear.
Definition 2. Quadratic Programming Problem
A Quadratic Programming Problem (or simply Quadratic Problem) is a linearly constrained mathematical optimization problem of a quadratic function. A QP is a particular type of NP. The quadratic function may be defined with respect to several variables, all of which may be subject to linear constraints. Considering a vector, a symmetric matrix , a real matrix , a real matrix , a vector and another vector , the goal of a QP is to determine the vector that minimizes a regular quadratic function of formsubject to constraints A_{ineq}x_{c} ≤ b_{ineq} and A_{eq}^{x}_{c} = b_{eq}.. The solution x_{c} to this kind of problem is found by many solvers seen in the literature, based on Interior Point algorithms, quadratic search, etc.
1.3. Notation
In this work, the set of nonnegative real number is denoted by , whist the set of nonnegative integers including zero is denoted by . The index set represents , with 0 ≤ a ≤ b. The identity matrix of size j is denoted as ; col {a, b, c} denotes the vectorization (collection) of the entries and diag{v} denotes the diagonal matrix generated with the line vector v.
The value of a given variable v(k) at time instant k + i, computed based on the information available at instant k, is denoted as v(k + i\k).
refers to the class of positive and strictly increasing scalar functions that pass through the origin. A given function is of class if f(0) = 0 and lim_{ξ→+∞}f(ξ) → +∞. A realvalued scalar function belongs to class if it belongs to class and it is radially unbounded (this is lim_{s→∞}ϕ(s) → +∞. A function belongs to class if, for each fixed and, for each fixed , β(s,·) is nonincreasing and holds for lim_{m→+∞}β(s,m) = 0.
denotes the set of all compact convex subsets of . A convex and compact set with nonempty interior, which contains the origin, is named a PCset. A subset of is denoted a polyhedron if it is an intersection of a finite number of half spaces. A polytope is defined as a compact polyhedron. A polytope can be analogously represented as the convex hull of a finite number of points in . A hyperbox is a convex polytope where all the ruling hyperplanes are parallel with respect to their axes.
Finally, consider two sets and . The Minkowski set addition is defined by A ⊕ B := {a+ba ∈ A, b ∈ B},while the Pontryag in set difference is defined by A ⊖ B := {a a ⊕ B ⊆ A}. The cartesian product between two sets is defined as .
2. PRELIMINARIES
In this Section, we detail how nonlinear processes can be described under a qLPV formalism; we also present some other formalities.
2.1. The Nonlinear System and its qLPV Embedding
We consider the following generic discretetime nonlinear system:
where represents the sampling instant, represents the system states, is the vector of control inputs and stands for the measured outputs of the process.
We begin by characterizing this process, which should satisfy the following key Assumptions:
Assumption 1.The admissible zone for the states is given by a 2norm upper bound on each entry x_{j}, this is:
Assumption 2.The admissible region for the control inputs is given by a 2norm upper bound on each entry u_{j}, this is:
Assumption 3.The nonlinear maps and are continuous and continuously differentiable with respect to x, i.e. class C^{∞}
Assumption 4.This nonlinear system is controllable in terms of x and y through the input trajectory u.
To represent any nonlinear system, as the one in Eq. (1), under a qLPV formalism, this last Assumption must be verified, since it is the LDI property that furnishes the settings for such representation.
The LDI property is as follows: suppose that, for each x, u and y and for every sampling instant k, there exists a matrix such that
where is the set within which the LDI property holds. Then, when there exists a matrix H(·) that verifies Eq. (4), the nonlinear model from Eq. (1) can be equivalently expressed as:
which is a qLPV formulation where represents the endogenous nonlinear function for the scheduling parameters. Note that ρ(k) is bounded and known online at each instant k, but generally unknown for any future instant k + j, .
We consider that the qLPV scheduling parameters have bounded rates of variations, this is: ρ(k + 1) = ρ(k) + ∂ρ(k + 1), being , ∀ k their variation rates. This is very reasonable for any practical application. To assume that ρ varies arbitrary implies in quite conservative control synthesis ^{[23]}.
3. PROBLEM SETUP
Regarding the qLPV embedded model in Eq. (5), we proceed by detailing how MPC can be applied to regulate and control this system. For the simplify of the reference tracking demonstrations, we drop the inputoutput energy transfer, i.e. D(ρ(k)) = 0. We note that the processes with D(ρ(k)) ≠ 0 can still be dealt with the proposed method, we no additional drawbacks.
The essential idea behind MPC is to consider a quadratic finitehorizon functional cost, which embeds the performance objectives of the system within this given horizon. The implementation resides in minimizing this cost with respect to a control signal sequence, using a model of the system in order to make predictions for the future variable values along the horizon. The optimization also includes the operational constraints of the process variables (admissibility region). Generically, we consider the following steadystate reference tracking performance cost:
where Q and R are positive definite weighting matrices and the pair (x_{r}, u_{r}) defines a known admissible steadystate reference target for the nonlinear system. The optimization cost J considers a prediction horizon of N_{p} steps and a positive terminal stage value V (x(k + Npk)) > 0.
The MPC framework considers a movingwindow strategy. Therefore, at each sampling instant k, since x(k), and ρ(k) are known, the corresponding optimization problem is solved, which gives the solution . This solution constitutes the following sequence of control inputs
whose first input u(kk) = I_{1}U_{k} is applied to the process. Then, the horizon slides forward and the procedure is updated. The complete optimization, at each sampling instant k, is given as follows:
where X_{f} and V(·) are the terminal ingredients, combined to ensure recursive feasibility of the algorithm (see Section Stability and Offline Preparations).
Due to Eq. (8), it follows that the future values for the qLPV scheduling variables ρ(k + i) are not known for any i ≥ 1. At each instant k, the optimization operates based on the knowledge of x(k) and u(k), which can be used to compute ρ(k) through the nonlinear qLPV proxy f_{ρ}(·). One could easily include this proxy into the optimization, making it also subject to ρ(k + i) = f_{ρ} (x(k + i), u(k + i)) together with the process model, but this would convert Eq. (8) into a Nonlinear Programming Problem, which is associated with numerical complexity issues (as previously discussed).
The NP execution is computationally unattractive^{[22]} because of this general nonlinear dependence of the predicted states on the control inputs and on previous states. Therefore, following the lines of previous works ^{[19,22]}, this paper pursues a fast implementation of the LPV MPC optimization procedure in Eq. (8), which means that we do not seek to analytically include the nonlinear qLPV scheduling proxy f_{ρ} (·) to the optimization, but rather to provide values for the complete evolution of the scheduling parameters along the prediction horizon, as if they were known (thus detaching the nonlinear dependency). This is, we aim to solve Eq. (8) based on x(k), ρ(k) and on the future ”scheduling sequence” vector , being
If the actual evolution of the scheduling sequence is as gives P_{k}, the MPC ensures perfect regulation. Furthermore, it is formulated as a Quadratic Programming Problem, which can be tackled for many timecritical applications with modern solvers. In fact, we approximate the NP solution by one which resides in a ”guess” for the scheduling sequence P_{k}, which attractively converges to the actual value of this vector as the procedure iterates. The solution to estimate P_{k} is based on a Moving Horizon Estimation algorithm, which is further detailed in Section The MHEMPC Mechanism.
We must proceed by providing some complementary Assumptions regarding this qLPV MPC optimization problem setup. For such, we denote as the evolution of the state values along the prediction horizon, this is:
Assumption 5.The qLPV scheduling proxy is setwise and vectorwise applicable, this is, it holds asand also as f_{ρ} (X_{k}, U_{k}). The first operation stand for the application of f_{ρ}(·) to the bounds of each entry set, while the later stands for the application of f_{ρ}(·) to each sample of the entry vectors.
Assumption 6.The application of the scheduling proxy to the admissible zone for the states and inputs is a subset of the scheduling set, this is:
Assumption 7.The admissible regionis a subset of the image of the inverse of the scheduling proxy domain, being f_{ρ}(·) bijective. This means that the inverse of the scheduling proxy always maps admissibility pairs (x, u) from admissible scheduling variables ρ. This is mathematically expressed as follows:
From the viewpoint of each sampling instant k, the scheduling sequence can be directly evaluated as:
where comprises the instantaneous states and the state evolution X_{k} until x(k + N_{p} – 1k):
which is directly given by with the last entry suppressed.
With the previous discussion in mind, we proceed by using the qLPV model from Eq. (5) and the definitions from Eqs. (7), (9) and (10) to analytically provide a solution to the state evolution which is explicitly dependent on the scheduling sequence.
For the LTI case, the state evolution X_{k}, departing from x(k), is expressed on a linear dependent basis w.r.t. x(k) and to the sequence of control inputs U_{k}, as follows:
Analogously, for the qLPV case, since linearity is retained through the inputoutput channels (i.e. from u to x), the state evolution can be given in a quite similar fashion, but with parameter dependence on P_{k} appearing on the transition matrices, this is:
where the parameter dependent matrices are given by^{1}:
In order to compute matrices and , some nonlinear operations should be performed. Anyhow, the procedure to compute them can be done completely outside the MPC optimization. In this form, the MPC receives, at each instant k, the following inputs: x(k), and (as well as the steadystate target given by x_{r} and u_{r}); then, solving the optimization problem in Eq. (8), it results in U_{k}, from which the first entry u(kk) is applied to the plant. For such goal, the MPC requires to internally explicitly minimize the cost function J(x, u, k) from Eq. (6), which can be written in the vector form as follows:
where H(P_{k}) is the Hessian of this cost function, g(·) its gradient and κ(·) an offset term. The notation Q̆ and R̆ denote the blockdiagonal version of these matrices, i.e. and , respectively, while and . The Hessian, gradient and offset terms are analytically given by:
3.1. Process constraints
The qLPV MPC problem solution proposed in this paper is formulated with respect to the scheduled state evolution equation, as gave Eq. (15), with H(P_{k}), g(·) and κ(·) being as passed as inputs to the resulting MPC optimization. This means that the MPC optimization does not treat state evolution X_{k} as optimization variables, but the whole problem is formulated singularly in terms of U_{k}.
For this reason, the admissibility process constraints and are conversely written in the following fashion, w.r.t. , , U_{k} and x(k), instead of u(k + j) and x(k + j):
The above formulations also facilitates the inclusion of slew rates on u, i.e. constraints on the control variation δu(k + i) = u(k + i) – u(k + i – 1), which can be done directly by adapting the vectorwise set .
The terminal constraint x(k + Npk) ∈ X_{f} is stated in terms of U_{k} as:
where
Additional output constraints are also easily formulated. If the process has some outputs y_{c} (which are not necessarily equal to y, but could be) that must be hard constrained, i.e. . Then, assume that these outputs can be described as:
In what follows, we take
Then, the additional constraint is formulated as follows^{2}:
thus
3.2. Reference tracking
Finally, before showing the proposed mechanism to guess P_{k} and solve the qLPV MPC problem, a comment must be made regarding reference tracking. The considered cost function J(x, u, k) from Eq. (6) (or its vector form of Eq. (16)) is set in order to minimize the variations from the desired setpoint target p_{r}= (x_{r}, u_{r}).
The majority of processes that require reference tracking, require it regarding the controlled outputs and not the states. This is, to ensure that y(k) tracks some steadystate value y_{r}. Since the controlled outputs in Eq. (5) are given by
we can find a linear (parameter varying) combination of the states x that, if tracked, ensures that y(k) → yr. We denote the output tracking target as , which is known.
Then, following the lines of previous reference tracking frameworks ^{[2426]}, we use an offline reference optimization selector, which is set to find the setpoint target p_{r} that abides to the constraints and ensures an output tracking of y_{r}. This nonlinear optimization procedure is as follows:
This procedure ensures some steadystate x_{r}= A(ρ_{r})x_{r}+ B(ρ_{r})u_{r} ) that abides to the states constraints and guarantees that the output tracking goal y_{r} is followed.
Note that this optimization procedure has a steadystate target point p_{r}= (x_{r}, u_{r}) as output, and not the full state and input trajectories towards this target.
The state reference selection problem can be solved online, at each sampling instant, if the output reference goal y_{r} changes over time. By doing so, an additional computational complexity appears, which can be smoothed if the scheduling parameter guess P_{k} is used instead of solving the nonlinear optimization itself. A full discussion on periodically changing reference tracking for nonlinear MPC has been recently presented ^{[27]}. The focus of this paper is constant reference signals, either given in terms of states x_{r} or outputs y_{r}
It is important to notice that, in order for the method to hold, the state reference x_{r} must be contained inside the terminal set of the MPC problem from Eq. (8). This ensures that the stability and recursive feasibility guarantees (as verified in Section Stability and Offline Preparations) hold.
4. THE MHEMPC MECHANISM
In the general qLPV embedding case of Eq. (5), the scheduling proxy f_{ρ}(·) is an arbitrary function of both state and input. For notation ease, we drop the control input dependency, using taking ρ(k) = f_{ρ}(x(k)). Anyhow, note that all that follows can be trivially extended to broader case.
The backbone idea of the method proposed in this paper follows the fashion of previous papers^{[19,22]}: to iteratively refine the predictions/guesses of the scheduling sequence P_{k} based on the (adjusted) state predictions . The main novelty of this paper is how the refining and estimation of P_{k} is done: in the prior, the scheduling sequence is taken, as each iteration, directly as gives Eq. (13), i.e. as a nonlinear operator upon a vector, which can be computationally difficult to track, according to the kind of nonlinearity present on the scheduling map f_{ρ}(·) ; in contrast, in this paper, P_{k} is taken according to a linear timevarying operator on , this is: . This linear operator derives from a MovingHorizon Estimation procedure, which proceeds by trying to match a fixedsize linear autoregressive model for P_{k}, being Θ_{k} the model parameters computed by the MHE at a given sampling instant k and P_{k1} the scheduling sequence guess at the previous instant.
A priori, the operation of this MHE has the computational complexity of a QP, which could be faster to evaluate than the Eqs. (13) with and , depending on the amount of nonlinearities present in f_{ρ}(·). Moreover, the proposed MHEMPC mechanism is able to provide convergence to the real scheduling sequence P_{k} faster than looping Eq. (13) to the MPC^{[22]}, as demonstrated through the experiment presented in Section Benchmark Example. The method follows:
Assumption 8.The scheduling map f_{ρ}(·) is algebraic.
Proposition 1.The scheduling map can be approximated by the following regression:
which can be given in compact form by:
beinga fairly easy map to compute with respect to numerical burden.
Proof. Indeed, due to Assumption 8, it is quite reasonable that Assumption 1 holds: any algebraic function of form ρ(k + 1) = f_{ρ}(x(k + 1)) can be Taylorexpanded to achieve a linear dependency on x with sufficiently small error, i.e.
This concludes the proof. □
The proposed procedure uses a MHE mechanism to estimate these parameter values a_{0,0} to a_{Np–1,Np–11}, at each sampling instant, concatenated as Θ_{k},through the following QP:
black
Essentially, this MHE scheme operates in order to find a parameter matrix Θ_{k} that makes the linear autoregressive equation yield the best match between the state evolution sequence and the qLPV model. Notice that the MHE algorithm only needs a few data from the previous step to find the parameter matrix Θ_{k}, being these the state predictions , the scheduling sequence P_{k–1} and the input vector U_{k}.
Figure 1 syntheses the proposed algorithm, which relies in a coordination between the MHE and the MPC optimization procedures. We must note that the MHE loop should operate until P_{k} converges to the actual value for the scheduling sequence, or until a certain stop criterion/heuristic threshold for the number of iterations is reached.
Proposed MHEMPC Scheme using a qLPV Model.
The proposed scheme is also detailed through the Algorithm below. Its application departs with an initial state evolution sequence X_{0}, that can be simply taken as constant/frozen evolution for the states and two initial scheduling sequences P_{0} and P_{–1}, which can be simply taken as if ρ(k) remained frozen along the whole prediction horizon, i.e. . The Algorithm also departs with a known terminal set condition Xf and a known target reference goal p_{r}. The Hessian, gradient and offset of J_{k} are taken as give Eqs. (17)(19), respectively.
4.1. Convergence property
In order to demonstrate the convergence of the proposed method (as in its implementation form given in Algorithm 1), we will proceed by verifying a wellknown results for Newton based SQPs from the literature ^{[2830]}. This is the same path followed in a previous paper ^{[22]}, which invokes established results to demonstrate that, under certain conditions, the MHEMPC mechanism that is solved at each iteration is equivalent to a quadratic subproblem used in standard Newton SQP. Therefore, a local convergence property is readily found.
Proposed qLPV MPC
Proposition 2.A quadratic subproblem program of SQP algorithms is derived by a secondorder approximation of the SQP optimization cost and a linearization of its constraints.
Proof. Found in ^{[28]}. □
For illustration purposes regarding this matter, we consider the following generic NP (as given in Definition 1):
This problem can be given as a quadratic subproblem directly, as follows:
where H_{fc} (x_{c}) denotes the Hessian of the optimization cost f_{c}(x_{c}) and ∇h_{j}(x_{c}) and ∇g_{i}(x_{c}) denote divergent operators. We note that this subproblem is evaluated at a given solution estimate x̅_{c} (at some given iteration), for which .
Regarding the proposed MHEMPC mechanism, we can easily show that if either simple Jacobian linearization or Linear Differential Inclusion are used to find a qLPV model (as in Eq. (5)) for the nonlinear system (as in Eq. (1)), then, the proposed mechanism iterates in equivalence to a Newton SQP subproblem. Notice how such subproblem in Eq. (27) is identical to the MHEMPC optimization given through the consecutive iterations of the MHE (Eq. (25)) and the MPC (Eq. (8)). The terminal constraint in the MPC optimization adds no convergence tradeoff.
Thence, it follows that if local convergence of the equivalent Newton SQP can be established, the proposed MHEMPC also yields convergence. The sufficient conditions for local convergence of a Newton SQP subproblem at x_{c} =x̅_{c}, as given by prior references ^{[2932]}, are that (i) the problem is set simply with equality constraints (not the MPC case) or that (ii) the subset of active inequality constraints are known before the optimization solution. The second condition is also not true for general MPC paradigms. However, one can iterate the subproblem until convergence is found at another point , as previously discussed ^{[22,33,34]}.
In practice, the proposed MHEMPC will not be set to freely iterate until the convergence of P_{k}. This is not desirable because the number of iterations needed for convergence may require more time than the available sampling period. Therefore, a stop criterion is added to the mechanism, so that iterations stop at a given threshold. A warmstart is also included by shifting the result regarding X_{k} and P_{k} from one sampling k as the initial guess for the optimization at k + 1, which ensures that the proposed algorithm reaches convergence after a few discretetime samples. We note that convergence of Newton SQP subproblems with warmstart have been assessed ^{[33,34]} and shown to be practicable for realtime schemes.
5. STABILITY AND OFFLINE PREPARATIONS
In this Section, we offer a Theorem to construct the terminal ingredients of the MPC algorithm: (a) the terminal set within which x(k + N_{p}k) is bounded to, and (b) the terminal offset cost V(x(k + N_{p}k)) minimized by the MPC.
Since the establishment of terminal ingredients toolkit ^{[2,35]} as the key way to ensure stability and recursive feasibility of statefeedback predictive control loops, MPC grew on both theory and industrial practice.
The usual approach with terminal ingredients resides in some ensuring that conditions are met by (a) the terminal set X_{f} and (b) the terminal cost V(x(k + N_{p}k)) with respect to a nominal statefeedback controller u(k) = k^{n}x(k), which is usually the unconstrained solution of the MPC problem. For the tracking case, the nominal feedback is given by u(k) = k^{n} (x(k) – x_{r}) (and so is the terminal constraint (x(k + N_{p}k) – x_{r}) ∈ X_{f} and the terminal cost V(x(k + N_{p}k) – x_{r})). Accordingly we develop a sufficient stability condition for the proposed MHEMPC mechanism in order to verify these conditions.
Firstly we consider that there exists a parameterdependent nominal statefeedback gain . For demonstration simplicity and notation lightness, we will proceed with x_{r} null^{3}.
This nominal controller is purely fictional, used to demonstrate stability and recursive feasibility properties of the proposed MHEMPC mechanism. Anyhow, it stands for the infinitehorizon LPV Linear Quadratic Regulator (LQR) solution, which verifies and the admissibility of x(k + i) and u(k + i – 1).
Of course, there is a complexity barrier to solve this problem, because the states have parametric nonlinearities that impact their trajectories (the qLPV scheduling parameters). Therefore, we determine this nominal feedback gain together with the terminal ingredients, which are also taken as parameterdependent on ρ. We consider, for regularity, an ellipsoidal set as the terminal constraint, which is given by:
This ellipsoid is centered at the origin and has a radius of α_{p}. Furthermore, this terminal set is a sublevel set of terminal cost V(·), which is taken as a Lyapunov function as follows:
This parameterdependent nominal feedback gain K^{n}(ρ)) and the parameter dependent terminal ingredients verbalized through the symmetric parameter dependent Lyapunov matrix P(ρ) are so that the following inputtostate stability Theorem is guaranteed.
Let Assumptions 4 and 7 hold. Assume that a nominal control law u = K^{n}(ρ)x exists. Consider that the MPC is in the framework of the optimization problem in Eq. (8), with a terminal state set given byX_{f}(ρ) and a terminal cost V(x,ρ). Then, inputtostate stability is ensured if the following conditions are hold:
(C1) The origin lies in the interior ofX_{f}(ρ);
(C2) Any consecutive state to x, given by (A(ρ) + B(ρ)K^{n}(ρ)) x lies withinX_{f}(ρ) (i.e. this is an invariant set);
(C3) The discrete algebraic Ricatti equation is verified within this invariant set, this is, ∀ x ∈ X_{f} (ρ(k)):
(C4) The image of the nominal feedback always lies within the admissible control input domain:.
(C5) The terminal setX_{f} (ρ) is a subset of the admissible state domain .
Assuming that the initial solution of the MPC problem , computed with respect to the initial state x(0), is feasible, the MPC algorithm is indeed recursively feasible, asymptotically stabilizing the state origin.
Proof. Provided in Appendix Proof of Theorem 1. □
In order to find some nominal statefeedback gain K^{n}(ρ), some terminal set X_{f} and some terminal offset cost V(·), an offline LMI problem is proposed in the sequel. This LMI problem is such that a P(ρ) positive definite parameterdependent matrix is found to ensure that the conditions of Theorem 1 are satisfied. Due to condition (C3), the LMI is solved over a sufficiently dense grid over ρ, consider its admissibility domain . We note that (C3) is a timevariant condition, which depends explicitly on ρ(k) and ρ(k + 1) due to the nature of the parameter dependent V(·).
This LMI problem is provided through the following Theorem, which aims to find the largest terminal set X_{f} that is invariant under the nominal control policy u(k) = K^{n}(ρ(k))x(k) for all k, while remaining admissible, i.e. , and Note that the largest ellipsoidal set as in the form of Eq. (28) is posed through the maximization of α_{p}.
Theorem 2.Terminal Ingredients^{[22,37]}
The conditions (C1)(C5) of Theorem 1 are satisfied if there exist a symmetric parameterdependent positive definite matrix a parameterdependent rectangular matrix and a scalar such that Y(ρ) = (P(ρ))^{–1}> 0, W(ρ) = K^{n} (ρ)Y (ρ) and that the following LMIs hold for alland, while minimizingα̂_{P}:
where I_{j} denotes the jthrow of the identity matrix. In LMI (31), it is given w.r.t. to an identity, while in LMI (32), it is given w.r.t. to an identity.
Proof. Refer to^{[37]}. □
We must note that the above proof demonstrates that the solution of the LMIs presented in Theorem 2 ensure a positive definite parameter dependent matrix P(ρ) which can be used to compute the MPC terminal ingredients V(·) and X_{f} such that inputtostate stability of the closedloop in guaranteed, verifying the conditions of Theorem 1. Furthermore, when the MPC is designed with these terminal ingredients, for whichever initial condition x(0) ∈ X_{f} it starts with, it remains recursively feasible for all consecutive discrete time instants k > 0.
Anyhow, Theorem 2 provides infinitedimensional LMIs, since they should hold for all and for all . To address this issue, one can handle the LMIs considering an sufficiently dense grid ^{[38]} of points in , for which the LMIs must be enforced. This solves the infinite dimension of the problem, which is converted into an n_{g}dimensional LMI problem, being n_{g} the number of grid points. For this solution to be practically implementable, continuity on matrices A(ρ) and B(ρ) should be verified. We must also note that parameterdependency on ρ may be dropped if the system is quadratically stabilizable, which is a conservative assumption.
6. BENCHMARK EXAMPLE
In this Section, we pursue the application of the proposed MHEMPC mechanism with terminal ingredients found through Theorem 2. For such, we consider the application of our control method upon a benchmark system, detailed in the sequel.
6.1. Continuouslystirred tank reactor
Consider the model of a ContinuouslyStirred Tank Reactor (CSTR) process, which consists of an irreversible, exothermic reaction, A → B, in a constant volume reactor cooled by a single coolant stream which can be modeled by the following equations:
where C_{A} is the concentration of A in the reactor, T is the temperature in the reactor, and T_{c} is the coolant temperature.
In this process, u = T_{c} is a control input, whereas C_{A} and T are measurable process variables. The considered model parameters and process constraints are reported in Table 1.
Model Parameters and Constraints
Parameter
Value/Set
Parameter
Value/Set
q
100 L min^{–1}
C_{Af}
1 mol L^{–1}
k_{0}
7.2 × 10^{10} min^{–1}
E/R
8750 K
H_{Δ}
–5 × 10^{5} cal mol^{–1}
ρC_{p}
239 cal L^{–1} K^{–1}
W
7 × 510^{4} cal min^{–1} K^{–1}
V
100 L
T_{f}
350K


C_{A}
∈ [0.03, 0.12] mol 1^{–1}
T
∈[440,460] K
Tc
∈ [200, 380]K


6.2. qLPV Embedding
Considering x = (C_{A}, T) as state variables, we obtain the following qLPV realization of the CSTR system from Eq. (33):
where ρ = f(x) denotes two scheduling variables, given as linear functions of each state, i.e. ρ_{1}= f_{1}(x_{1}) and ρ_{2}=f_{2}(x_{2}). This continuoustime model is Eulerdiscretized using a sampling period of T_{s}= 30ms.
6.3. Control goal and tuning
Considering an arbitrary initial condition x_{0} given within the state admissibility set , the proposed controller is set to steer the state trajectories to a known state reference x_{r}. For such, we use identity weights in the MPC cost J(x, u, k). Complementary, we use a prediction horizon of N_{p}= 8 steps. This prediction horizon was chosen in accordance with prior literature using the same nonlinear CSTR benchmark^{[39]}.
The MHE scheme, which is used to estimate the future scheduling behaviour P_{k} at each sampling instant k, is set to operate with a threshold loop barrier of 3 loops, which is a verified sufficient bound to induce convergence (refer to the discussion in Section Convergence Property).
6.4. Simulation results
Considering a realistic nonlinear CSTR model, we obtain simulation results to demonstrate the effectiveness of the proposed control scheme. The following results were obtained in a 2.4 GHz, 8 GB RAM Macintosh computer, using Matlab, yalmip and Gurobi solver.
First, we show how the MHE operation is able to accurately predict the behaviour of the scheduling trajectory, as depicts Figure 2. At each instant k, the MHE provides estimation for P_{k}, which is composed of the following N_{p} entries of the scheduling variables. In this Figure, we observe the real scheduling variables ρ_{1} and ρ_{2} (dotdash black line) and the estimates P_{k} provided at different samples (coloured xmarked lines). Within some samples, we can see that the predicted trajectory converges to the real one, which confirms the effectiveness of the MHE operation.
MHE: scheduling trajectory estimates.
Based on the scheduling behaviours predicted by the MHE loop, the predictive controller determines the control input (Figure 4) in order to drive the system states from x_{0} to the reference goal x_{r}. The corresponding state trajectories are depicted in Figure 3, which also shows the state admissibility set and the terminal set constraint X_{f} (a parameterdependent ellipsoid generated via Theorem 2). As one can see, the behaviour of the process variables is a smooth trajectory towards x_{r}.
State trajectories, state admissibility set, terminal set.
MPC: control input.
Finally, we demonstrate the dissipating properties of the proposed control scheme. In Figure 5, we show the evolution of MPC stage cost J over time. As expected, J decays and converges to the origin, which verifies the dissipation properties required by Theorem 1.
MPC: stage cost.
7. CONCLUSIONS
In this paper, a new method for the fast, realtime implementation of Nonlinear Model Predictive Control is proposed. The method provides a nearoptimal, approximated solution, which is found through the online operation of sequential Quadratic Programming Problems. The main necessary argument to develop the method is that the nonlinear process should be described by quasiLinear Parameter Varying model, for which the embedding is ensured through a scheduling proxy. Then, the online operations resides in the consecutive operation of the MPC program together with a MovingHorizon Estimation scheme, which is used to match the future values of the scheduling proxy along the prediction horizon, which are unknown. Inputtostate stability and recursive feasibility properties of the algorithm are ensured by parameterdependent terminal ingredients, which are computed offline. Using a benchmark example, the method is tested. We highlight that it proves itself more effective for stronger nonlinearities in the qLPV scheduling proxy, for which the MHE scheme operates faster than the application of the scheduling proxy upon each entry of the future state variables, as in many other techniques. For future works, the Authors plan on assessing the issue of periodicallychanging (possibly unreachable) output reference signals.
Proof of Theorem 1
Consider an initial state condition x(k_{0}) with a scheduling parameter p(k_{0}) (transition map). Assume that this initial condition generates a feasible optimal control sequence
The next feasible sequence, which might be possibly suboptimal, at instant k_{0} + 1, is denoted:
Notice that the last entry of this second sequence u(k_{0}+ N_{p}k_{0} + 1) is equal to the feedback of the terminal state from the first sequence, this is:
Then, let us^{4} evaluate the evolution the MPC cost function J(·) and its decay between time instants k_{0} and K_{0} + 1:
Assuming that was applied to the plant at instant k_{0} (this input is the first entry of the optimal solution at time K_{0}), it follows that:
From (C3), we can pursue with the negativeness of the following term:
Thus, substituting this inequality into Eq. (35) yields:
Since Q and R are positive defined matrices by construction (they are the MPC tuning weights) and and x(k_{0} + 1) are known^{5}, it follows that:
which means that the MPC cost function decays along k.
From the optimality of the solution , it follows that the cost constructed with this sequence holds as a lowerbound with respect to J(k_{0} + 1), this is . Then, using this argument on ΔJ(k_{0}), we arrive at:
which proves that is a Lyapunov function and x will converge to the origin (due to (C1)), as long as the initial condition provides, in fact, a feasible starting point. We note that (C2) is necessary to map a feasible x(k_{0}+ 1k_{0}). Conditions (C4) and (C5) are necessary to ensure admissible trajectories regarding x and u.
DECLARATIONSAuthors’ contributions
All authors contributed equally.
Availability of data and materials
Data will be made available upon email request to the corresponding author.
Financial support and sponsorship
M. M. Morato is partially supported by CNPq project 304032/2019 – 0. V. Stojanovic thanks the Serbian Ministry of Education, Science and Technological Development for support (grant No. 451039/202114/200108).
Conflict of interest
Both Authors declare no potential conflict of interests.
^{2}Note that the last y_{c}(k + N_{p}k) is not constrained due to unavailability of ρ(k + N_{p}) inside P_{k}. This issue is amendable by taking a longer scheduling prediction guess P_{k}, but out of the scope.
^{3}The tracking equivalency is easily done with and by computing the qLPV model with the states dynamics are given with respect to x_{tracking}(k) = x(k) – x_{r}.
^{4}We use compact notation for brevity. Herein, we use ρ := ρ(k_{0}+ 1).
^{5}Note that .
CamachoEFBordonsCMayneDQRawlingsJBRaoCVScokaertPOMConstrained model predictive control: Stability and optimality.AllgöwerFZhengACamachoEFBordonsCGrosSZanonMQuirynenRBemporadADiehlMFrom linear to nonlinear MPC: bridging the gap via the realtime iteration.QuirynenRVukovMZanonMDiehlMAutogenerating microsecond solvers for nonlinear MPC: a tutorial using ACADO integrators.EnglertTVölzAMesmerFRheinSGraichenKA software framework for embedded nonlinear model predictive control using a gradientbased augmented lagrangian approach (GRAMPC).OhyamaSDateHRathaiKMMSenameOAlamirMGPUbased parameterized nmpc scheme for control of half car vehicle with semiactive suspension system.HoffmannCWernerHA survey of linear parametervarying control applications validated by experiments or high fidelity simulations.MoratoMMNormeyRicoJESenameOModel predictive control design for linear parameter varying systems: A survey.BoydSEl GhaouiLFeronEBalakrishnanVAbbasHSTothRPetreczkyMMeskinNMohammadpourJEmbedding of nonlinear systems in a linear parametervarying representation.KunzKHuckSMSummersTHCisnerosPablo SGWernerHerbertWide range stabilization of a pendubot using quasiLPV predictive control.AlcaláEPuigVQuevedoJLPVMPC control for autonomous vehicles.MateSKodamanaHBhartiyaSNatarajPSVA stabilizing suboptimal model predictive control for quasilinear parameter varying systems.MoratoMMNormeyRicoJESenameONovel qLPV MPC design with leastsquares scheduling prediction.MoratoMMNormeyRicoJESenameOSuboptimal recursively feasible linear parametervarying predictive algorithm for semiactive suspension control.CisnerosPSGVossSWernerHCisnerosPGWernerHFast nonlinear MPC for reference tracking subject to nonlinear constraints via quasiLPV representations.CisnerosPSWernerHNonlinear model predictive control for models in quasilinear parameter varying form.JungersMCaunRPOliveiraRCLFPeresPLDModel predictive control for linear parameter varying systems using pathdependent lyapunov functions.LimonDFerramoscaAAlvaradoIAlamoTCamachoEFLimonDFerramoscaAAlvaradoIAlamoTNonlinear MPC for tracking piecewise constant reference signals.MorariMMaederUNonlinear offsetfree model predictive control.KöhlerJMüllerMAAllgöwerFA nonlinear tracking model predictive control scheme for dynamic target signals.QiLSuperlinearly convergent approximate newton methods for lc^{1} optimization problems.WeiZLiuLYaoSThe superlinear convergence of a new quasinewtonsqp method for constrained optimization.IzmailovAFSolodovMVOn attraction of linearly constrained lagrangian methods and of stabilized and quasinewton sqp methods to critical multipliers.BoggsPTTolleJWKearsleyAJBoggsPTTolleJWSequential quadratic programming for largescale nonlinear optimization.DiehlMBockHGSchlöderJPA realtime iteration scheme for nonlinear optimization in optimal feedback control.HouskaBFerreauHJDiehlMAn autogenerated realtime iteration algorithm for nonlinear MPC in the microsecond range.MichalskaHMayneDQRobust receding horizon control of constrained nonlinear systems.DuanGRYuHHMoratoMMNormeyRicoJSenameOShortsighted robust lpv model predictive control: Application to semiactive suspension systems.WuFA generalized LPV system analysis and control synthesis framework.ChenHKremlingAAllgöwerFNonlinear predictive control of a benchmark CSTR.