Marcelo Menezes Morato
email: marcelomnzm@gmail.com

Victor M Cunha

Tito L M Santos

Julio E Normey-Rico

Olivier Sename

A Robust Nonlinear Tracking MPC using qLPV Embedding and Zonotopic Uncertainty Propagation

Keywords: Tracking Model Predictive Control, Linear Parameter Varying Systems, Zonotopes

published or not. The documents may come

Introduction

Context

Model Predictive Control (MPC) is indisputably well established. Nevertheless, the application of robust Nonlinear MPC (NMPC) is not trivial and comes with increased numerical burden, which may be an impediment for real-time applications. The majority of stabilizing NMPC schemes ensure regulation of the closed-loop dynamics to a fixed target [START_REF] Boccia | Stability and feasibility of state constrained MPC without stabilizing terminal constraints[END_REF]. Accordingly, asymptotic stability and constraints satisfaction are usually guaranteed with terminal ingredients, which verify invariance conditions in the neighborhood of the operation target. Nevertheless, this design method is not valid for set-point changes, since feasibility may be lost [START_REF] Köhler | A nonlinear model predictive control framework using reference generic terminal ingredients[END_REF].

Therefore, there has been an increasing focus on NMPC schemes for Tracking, considering time-varying setpoints. We highlight the following concrete Tracking NMPC tools for piece-wise constant reference signals: (i ) the use of artificial reference signals, as proposed by [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF], which allows less conservative terminal constraints and ensures feasibility is not lost; and (ii ) the generalisation of the prior for the case of periodic reference signals with the use of terminal equality constraints or online optimized terminal sets, as proposed by [START_REF] Köhler | A nonlinear tracking model predictive control scheme for dynamic target signals[END_REF].

Motivations

In this paper, we also address the issue of tracking possibly unreachable output target signals using statefeedback NMPC, closely building upon these previous papers. We address the following disadvantages of the prior methods: (i ) they require the online solution of Nonlinear Programs (NPs), which are numerically ex-pensive and not viable for time-critical applications 1 ; and (ii ) no prediction model mismatches are considered, which should be included in any real application.

We must mention that, in parallel to the theoretical establishment of NMPC, the Linear Parameter Varying (LPV) toolkit has been brought to focus [START_REF] Sename | Robust control and linear parameter varying approaches: application to vehicle dynamics[END_REF]. LPV models can be used to represent nonlinear dynamics with the use of known, bounded scheduling parameters ρ [START_REF] Mohammadpour | Control of linear parameter varying systems with applications[END_REF]. Recent advances on NMPC algorithms have been presented by exploiting the quasi-LPV (qLPV) embedding, thus enabling fast implementation, e.g. [START_REF] Morato | Model predictive control design for linear parameter varying systems: A survey[END_REF] and references therein. The elegance of the qLPV approach is that the "full-blown" nonlinear predictions are replaced by linear laws.

Contributions

With respect to this context, we develop a NMPC algorithm for tracking, using nominal predictions obtained through qLPV embedding. Up to our best knowledge, qLPV-embedding NMPC has only been formalized for regulation purposes [START_REF] Menezes Morato | Novel qLPV MPC design with least-squares scheduling prediction[END_REF][START_REF] Cisneros | Nonlinear model predictive control for models in quasi-linear parameter varying form[END_REF].

Accordingly, as suggested in [START_REF] Köhler | A nonlinear tracking model predictive control scheme for dynamic target signals[END_REF], the proposed solution incorporates robustness qualities against model uncertainties and additive load disturbances, using constraint tightening under the framework of [START_REF] Tito | A constraint-tightening approach to nonlinear model predictive control with chance constraints for stochastic systems[END_REF]. Based on the one-step-ahead disturbance propagation, we enforce the satisfaction of the performance requirements by binding the prediction error within zonotope extensions [START_REF] Victor | Nonlinear robust predictive control with distrubance propagation via zonotopes[END_REF].

The main novelties of this work are:

(1) A method for the recursive extrapolation of scheduling parameters is developed (Sec. 3). The extrapolation serves to compute the qLPV model at each sampling instant, with bounded prediction errors. Complementary, these bounds are used to build the zonotopes that bind the uncertainty propagation along the prediction horizon2 . (2) We also offer robust parameter-dependent terminal ingredients, which ensure recursive feasibility and stability (Sec. 4). A complementary optimisation is proposed for the choice of the artifical reference with relivied complexity. (3) Finally, using a benchmark example, we demonstrate that the numerical complexity of the proposed online algorithm is, on average, comparable to Quadratic Programs (QPs), being much faster than the NPs from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF][START_REF] Köhler | A nonlinear tracking model predictive control scheme for dynamic target signals[END_REF] (Sec. 5). This property is made viable due to the qLPV model realisation, which enables linear predictions at each sampling period.

Preliminary discussions and the overall problem setup is presented in Sec. 2. General conclusions are drawn in Sec. [START_REF] Köhler | A nonlinear model predictive control framework using reference generic terminal ingredients[END_REF].

Remark 1 With respect to the Authors' previous works, we note that this paper enhances and generalises the robust regulation algorithm from [START_REF] Menezes Morato | Robust nonlinear predictive control through qlpv embedding and zonotope uncertainty propagation[END_REF] for the case of time-varying reference signals. The recursive extrapolation procedure is a novel adaptation of the method from [START_REF] Menezes Morato | Novel qLPV MPC design with least-squares scheduling prediction[END_REF], now with bounded residual.en

Notation. The index set N [a,b] represents {i ∈ N | a ≤ i ≤ b}, with 0 ≤ a ≤ b.
The identity matrix of size j is denoted as I j ; I j,{i} denotes the i-th row of I j ; col{•} denotes the vectorisation of the entries and diag{v} denotes the diagonal matrix generated with the line vector v. 1 n×m stands for the n × m vector of unit entries. 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). K refers to the class of positive and strictly increasing scalar functions that pass through the origin. A C 1 function f is such that it is differentiable with continuous derivatives. In this case, ∇ T f : R m → R n×m denotes its Jacobian matrix. Consider sets A, B ⊂ R n , C ⊂ R m and a matrix R ∈ R n×m . The Minkowski set addition is defined by

A ⊕ B := {a+b | a ∈ A , b ∈ B}, while the Pontryagin set difference is defined by A B := {a | a ⊕ B ⊆ A}. A linear mapping is RA = {y ∈ R n : y = Ra, a ∈ A}, while the Cartesian product holds as A × C = {z ∈ R n+m : z = (a T c T ) T , a ∈ A, c ∈ C}.The unitary m- dimensional box is denoted B m ∞ = {ξ ∈ R m : ||ξ|| ∞ ≤ 1}.
The set of real compact intervals is given by I = {[a, b], a, b ∈ R, a ≤ b}. An interval matrix J ∈ I n×m has mid(J) and rad(J) denoting its middle point and radius, respectively. • denotes the 2-norm.

Problem Statement

System Description

We consider the following qLPV-embedding of a nonlinear system:

x(k + 1) = A(ρ(k))x(k) + B(ρ(k))u(k) + w(k) y(k) = C(ρ(k))x(k) + D(ρ(k))u(k). (1)
where x ∈ R nx represent the states, u ∈ R nu the inputs, and y ∈ R ny the outputs. We consider that the additive disturbance w ∈ R nx is bounded to a compact set with the origin at its interior, in such a way that

w k ∈ W ⊆ R nx , w ≤ w.
Remark 2 We note that any nonlinear system can be embedded to a qLPV realisation, as long as (exact, linear, convex, convex-concave) differential inclusion is satisfied. This is detailed in [START_REF] Shamma | An overview of LPV systems[END_REF][START_REF] Morato | Model predictive control design for linear parameter varying systems: A survey[END_REF] and experimentally exemplified in [START_REF] Hoffmann | A survey of linear parameter-varying control applications validated by experiments or high-fidelity simulations[END_REF].

For set-operation simplicity, hard compact polyhedral constraints on state and input (x(k), u(k)) ∈ Z = X × U are considered, which define admissible operation of system given by:

X := {x ∈ R nx : H x x ≤ h x } U := {u ∈ R nu : H u u ≤ h u }. (2)
Outer box-type constraints3 are implied:

|x j | ≤ x j , ∀j ∈ N [1,nx] and |u j | ≤ u j , ∀j ∈ N [1,nu] . Y := {y ∈ R ny | H y y ≤ h y } = h(Z)
is a compact and convex set which defines the possible outputs y mapped by h(x, u

) = A(f ρ (x))x + B(f ρ (x))u with (x, u) ∈ Z.
This qLPV model is scheduled by the state-dependent time-varying parameter ρ(k) = f ρ (x(k)) ∈ P ⊆ R nρ , which is bounded and measured online (known). The nonlinear scheduling map f ρ (x) is algebraic, class C 1 for all x ∈ X , and set-and vector-wise applicable 4 .

Assumption 1

The nonlinear scheduling parameter map f ρ : X → P agrees to a local Lipschitz condition around any arbitrary point x ∈ X , this is:

f ρ (x) -f ρ (x) ∞ ≤ γ ρ (x -x) ∞ , ∀ x, x ∈ X , (3)
where the smallest constant γ ρ that satisfies Eq. (3) is known as the Lipschitz constant for f ρ (•).

We consider

P := {ρ ∈ R nρ : ρ j ≤ ρ j ≤ ρ j , j ∈ N [1,nρ] }.
We stress that the future scheduling parameters ρ(k + j), ∀j ∈ N [1,∞] are unknown at time instant k. Therefore, for the application of MPC, estimates ρ(k + j|k) are used, as further detailed in Sec. 3.

Assumption 2

The scheduling variables ρ(k) exhibit a bounded rate of variation. This is: δρ(k + 1) = (ρ(k + 1) -ρ(k)) ∈ δP, with:

δP := δρ j ∈ R : δρ j ≤ δρ j ≤ δρ j , ∀j ∈ N [1,np] . ( 4
)
Remark 3 Bounds on scheduling parameters' variations δP can be obtained based on their dependence on the states ρ(k) = f ρ (x(k)). Specifically, this can be done by placing bounds

5 on δρ = f ρ (A(f ρ (x))x+B(f ρ (x))u+ w)-f ρ (x) for (x, u), (A(f ρ (x))x+B(f ρ (x))u+w, u) ∈ Z
and w ∈ W, either by interval arithmetics or optimisation.

Closed-Loop Paradigm

Through the sequel, we assume that the states x(k) are measurable at each sampling instant. Therefore, in order to limit the disturbance propagation, we consider the following LPV state-feedback control law:

u(k) = v(k) + K(ρ(k))x(k), (5)
where the virtual input v(k) is used for the MPC synthesis and the parameter-dependent feedback gain K(ρ(k)) is chosen in order to attenuate the propagation of disturbances. We assume that this gain K(ρ(k)) is structurally known 6 . The closed-loop dynamics are thus given by:

x(k + 1) = (A(ρ) + B(ρ)K(ρ))x(k) + B(ρ)v(k) + w(k) = A π (ρ)x(k) + B(ρ)v(k) + w(k) , (6)
and the process constraints (x(k), u(k)) ∈ Z can be expressed in terms of state and virtual input by (x(k), v(k)) ∈ Z π , where:

Z π := z ∈ R n+nu : H x 0 H u K(ρ) H u z ≤ h x h u . (7)
Given an initial condition x(k) ∈ R nx , the closed-loop trajectories for the following j steps through the qLPV embedding are given by exploiting Eq. ( 6), which gives:

x(k + j) := φ j (x(k), v [k+j-1] , w [k+j-1] , ρ [k+j-1] ), (8)
being

v [k+j-1] , w [k+j-1]
, and ρ [k+j-1] the vectors of future inputs, disturbances and scheduling parameters from sample k to k + j -1, respectively. Note that this is a closed-loop model due to the stabilising feedback K(ρ)x, with matrix A π (ρ) structurally known, being v an external input (the MPC law).

Both the future disturbances (due to stochastic nature) and future scheduling parameters (due to future statedependence) are unknown from the MPC viewpoint at the sampling instant k. Therefore, we consider the following nominal prediction model:

x(k + j|k) = φ j (x(k), v [k+j-1] , 0, ρk [k+j-1] ), (9)
where the future disturbances are presumably null and the scheduling parameters estimates ρk [k+j-1] =

6 Consider an affine gain K(ρ) = nρ i=1 Kiρi. In this case, Ki are assumed to be known. col( r ho(k|k), . . . , r ho(k + j -1|k)) are provided by the extrapolation algorithm from Sec. 3 at sampling instant k. This model is equivalent to using:

x(k + 1|k) = A π ( r ho(k))x(k) + B(ρ)v(k|k), x(k + 1|k) = A π (ρ(k + 1|k))x(k + 1|k) + B(ρ)v(k + 1|k),
and so forth, with estimates ρ generated by the extrapolation algorithm detailed in Sec. 3.

Disturbance Propagation

Given a prediction horizon N p , we consider a group of compact sets E(j), j ∈ N [0,Np-1] that bounds the difference between the scheduling parameter predictions made in k and k+1. That is, we have ρ(k+j|k+1)-ρ(k+j|k) ∈ E(j), ∀j ∈ N [0,Np-1] 7 .

Due to the mismatch between nominal predictions in Eq. ( 9) and real system trajectories in Eq. ( 8), we consider one-step-ahead disturbance propagation sets S(j), j ∈ N [0,Np] [START_REF] Tito | A constraint-tightening approach to nonlinear model predictive control with chance constraints for stochastic systems[END_REF] in order to guarantee recursive feasibility of the MPC strategy.

Definition 1 One-step-ahead disturbance propagation sets S(j), j ∈ N [0,Np] are compact sets that satisfy: (1) S(0) contains W, and (2) For all

x a , x b ∈ R n , v ∈ R nu , ρ a , ρ b ∈ P and j ∈ N [1,Np] , with (x a , v) ∈ Z π (S(j-1)× {0}), (x b -x a ) ∈ S(j-1), and (ρ b -ρ a ) ∈ E(j-1), it fol- lows that (A π (ρ b )x b + B(ρ b )v)-(A π (ρ a )x a + B(ρ a )v) ∈ S(j).
Taking x a = x(k + 1|k) and x b = x(k + 1), we have, by induction from Def. 1, that x(k + j|k + 1) ∈ x(k + j|k) ⊕ S(j -1), ∀j ∈ N [START_REF] Boccia | Stability and feasibility of state constrained MPC without stabilizing terminal constraints[END_REF]Np+1] for any admissible sequence of inputs v [k+j-1] and scheduling parameter predictions ρk [k+j-1] and ρk+1

[k+j] .

The sets S(j) thus bound the difference between the predictions made in k and k + 1, and therefore can be used to guarantee recursive feasibility and constraint satisfaction of the MPC based on nominal predictions.

In order to compute zonotopic reachable sets S z (j), ∀j ∈ N [0,Np] that satisfy Def. 1, we develop Theorem 1, based on Lemma 1.

Lemma 1 Consider a centered zonotope X = M B ng ∞ ⊆ R m , an interval matrix J ∈ I n×m , and a zonotope family Z = JX = {Jx, J ∈ J, x ∈ X}. A zonotopic inclusion is defined by:

(Z) := mid(J)X ⊕ P B ng ∞ , (10)
7 Being Np the prediction horizon, at time instant k, we only need to predict ρ(k + j|k), j ∈ N [0,Np-1] in order to obtain

x(k + j|k) for j ∈ [1, Np]. We thus consider ρ(k + Np|k) = ρ(k + Np -1|k), ∀k ≥ 0.
where P is a diagonal matrix such that:

P ii = ng j=1 m k=1 rad(J) ik |M kj |, ∀i ∈ N [1,n] . (11)
It holds that Z ⊆ (Z).

Proof 1 Refer to [START_REF] Brenner S Rego | Guaranteed methods based on constrained zonotopes for set-valued state estimation of nonlinear discrete-time systems[END_REF].

Theorem 1 Consider Z ⊆ R nx+nu , S 0 ⊆ R nx zono- topes and A, ∆ A (j) ∈ I nx×nx , ∆ B (j) ∈ I nx×nu inter- val matrices satisfying W ⊆ S 0 , Z π ⊆ Z, A π (ρ) ∈ A, A π (ρ b (j)) -A π (ρ a (j)) ∈ ∆ A (j) and B(ρ b (j)) - B(ρ a (j)) ∈ ∆ B (j) for all ρ, ρ a (j), ρ b (j) ∈ P and j ∈ N [1,Np-1]
, with ρ b (j) -ρ a (j) ∈ E(j -1). The zonotopes S z (j), j ∈ N [0,Np] defined recursively by S z (0) = S 0 and

S z (j) = V(j) ⊕ (AS z (j -1)), j ∈ N [1,Np] , (12)
where V(j) is given by V

(j) = ∆ A (j) ∆ B (j) Z π , satisfy Def. 1. Proof 2 The first condition is satisfied by design. Con- sidering x a , x b ∈ R nx , v ∈ R nu , ρ a , ρ b ∈ R nρ
as given by the second condition for some j ∈ N [1,N p] , and

∆ j = (A π (ρ b )x b + B(ρ b )v) -(A π (ρ a )x a + B(ρ a )v),
we have:

∆ j = (A π (ρ b ) -A π (ρ a ))x a + A π (ρ b )(x b -x a ) + (B(ρ b ) -B(ρ a ))v ∈ ∆ A (j) ∆ B (j) Z π ⊕ AS z (j -1)
⊆ V(j) ⊕ (AS z (j -1)) = S z (j).

Therefore, the sets S(j) satisfy Def. 1.

Remark 4

In the case of qLPV systems with A π (ρ) and B(ρ) affine8 on ρ, it follows that A π (ρ b (j)) -A π (ρ a (j)) = A π (ρ b (j) -ρ a (j)) and B(ρ b (j)) -B(ρ a (j)) = B(ρ b (j) -ρ a (j)), with A π (•) and B(•) being linear mappings. Then, the interval matrices ∆ A (j) and ∆ B (j) can be computed directly from E(j). In the case of non-affine models, interval arithmetic or optimisation can be used to obtain the interval matrices ∆ A (j) and ∆ B (j) from P and E(j).

Remark 5 Due to the zonotope inclusion and Minkowsky sum, the number of generators of the zonotopes S(j) increases for each iteration. Methods for complexity reduction can be used to restrict the number of generators of each S(j) to a predefined value [START_REF] Joseph K Scott | Constrained zonotopes: A new tool for set-based estimation and fault detection[END_REF].

Remark 6

We stress that the qLPV model in Eq. ( 6) requires a linear nominal prediction model, as in Eq. ( 9), which changes at each sampling instant k, since it is based on an estimation of the future scheduling parameters. Therefore, the propagation of disturbances along the horizon is crucial, but may yield conservative sets S(j) due to the changes in ρk [k+j -1] . Nonlinear predictions can result in less conservative sets S(j), but also result in a higher online computational burden.

Remark 7 A trivial choice for ρk [k+j-1] would be ρk [k+j-1] = (ρ(k), . . . , ρ(k))
, amounting to a frozen estimation of the scheduling parameters. This would result in E(j) = δP, simplifying calculations. Nonetheless, the application of extrapolation algorithms, as the one described in Sec. 3, allow for less variation between ρk

[k+j-1] and ρk+1

[k+j] and thus smaller sets E(j), reducing the conservatism of the disturbance propagation.

Admissible Equilibrium Points

In this paper, we consider a state-feedback MPC strategy in order to steer the system states to a steady-state condition x r which implies in an output reference tracking goal y r . Therefore, we use the MPC algorithm to minimize the variations of (x, u) from the desired set-point targets (x r , u r ), which implies in the output target y r .

Following the lines of [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF], we consider that there exists a linear (parameter varying) combination of the states x and inputs u that ensures y(k) → y r . Through the sequel, we assume that for all y r ∈ Y there exists a pair z r = (x T r , u T r ) T which satisfies

(I -A(f ρ (z r ))) -B(f ρ (z r )) C(f ρ (z r )) D(f ρ (z r )) z r = 0 nx y r . ( 13
)
We consider the set of admissible steady-state outputs Y a ⊂ Y as the set of outputs generated by feasible states and inputs that satisfy the equilibrium condition, i.e., Eq. ( 13) is satisfied for some z r ∈ Z, with f ρ (z r ) ∈ P. More specifically, we have

Y a := {y ∈ Y | y = [C(f ρ (z r )) D(f ρ (z r ))]z r , z r ∈ Z Z ξ , z r = [A(f ρ (z r )) B(f ρ (z r ))]z r },
with Z ξ being a auxiliary nonempty set with the origin in its interior, included so that the frontier of Z is excluded (Z Z ξ ⊂ Z).

We assume that Eq. ( 13) implies in an unique correspondence between the steady-state output goal y r and the pair z r . Thus, there exist locally Lipschitz continuous functions g x : Y → R nx and g u : Y → R nu such that x r = g x (y r ) and u r = g u (y r ) for all y r ∈ Y a . Complementary, we assume that there exists a continuous functions g ρ : Y → R nρ that maps the equilibrium scheduling parameter ρ r = g ρ (y r ). Note that g ρ (y r ) := f ρ (g x (y r ), g u (y r )) and g v (y r ) = g u (y r ) -K(g ρ (y r ))g x (y r ). For such unique correspondence to exist, the left hand-side matrix in Eq. ( 13) must be non-singular for all z r ∈ Z [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF], which means there must be no redundant inputs or outputs as well as no integrator nodes in the qLPV dynamics in Eq. (1).

Tracking NMPC Design

Consider admissible time-varying piecewise constant output references y r ∈ Y a . Then, the nonlinear process must be controlled in such way that the state trajectories are suitably steered to x r , with a steady-state control input u r , despite disturbances and the discrepancy between the qLPV prediction model ( 9) and the real system trajectories [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF]. We consider a prediction horizon of N p steps.

As discussed in the prequel, we apply a robust NMPC scheme, for which contracted constraints are used. Considering an initial constraint set Z π (0) = Z π , the following sets for j ∈ N [START_REF] Boccia | Stability and feasibility of state constrained MPC without stabilizing terminal constraints[END_REF]Np] are iteratively taken as:

Z π (j + 1) = Z π (j) (S(j) × {0}) . (14)
Therefore, at each sampling instant k, we measure the state x(k), compute the scheduling parameter ρ(k), estimate the scheduling sequence ρk [k+Np-1] , and solve the following optimisation problem, which embeds the performance objectives of the system, as well as the operational constraints:

min v [k,k+Np-1] Np-1 j=0 (x(k + j|k) -x r , u(k + j|k) -u r )) +V (x(k + N p |k) -x r ) , (15)
s.t. :              x(k + j + 1|k) = A π (ρ(k + j|k))x(k + j|k) +B(ρ(k + j|k))v(k + j|k), j ∈ N [0,Np-1] , (x(k + j|k), v(k + j|k)) ∈ Z π (j), j ∈ N [0,Np-1] , x(k + N p |k) ∈ X f , . where (•, •) is a quadratic stage cost, V (•) is a terminal cost, x r = g x (y r ), u r = g u (y r
) and X f is a terminal robust invariant set (more discussions on this matter are presented in Sec. 4). Note that u(k + j|k) is implictly defined through Eq. ( 5). From the optimal solution of this problem v [k,k+Np-1] , we apply the first entry v (k|k) to the process according to Eq. ( 5).

A Recursive Extrapolation Algorithm

In order to estimate the scheduling sequence ρ[k+Np-1] , we use a rather simple extrapolation mechanism. The proposed method is an adaption of the framework in [START_REF] Menezes Morato | Novel qLPV MPC design with least-squares scheduling prediction[END_REF], now using a first-order Taylor expansion of the scheduling proxy f ρ (x(k)) around the state deviation instead of the Least-Squares rule. The main novelty is that the generated residual is bounded.

We denote δx(k + j) = x(k + j + 1) -x(k + j) as the incremental state deviation. This variable is bounded 9 to a compact and convex box-type set

δX := δx ∈ R nx : |δx j | ≤ δx j , ∀j ∈ N [1,nx] .
The method is as follows: consider that the static scheduling map f ρ (x) can be approximated by:

f ρ (x) = f ρ (x)| x + ∂f ρ (x) ∂x x (x -x) + ξ ρ , (16)
being x the expansion point and ξ ρ a residual noise which inherits the discrepancy between the real static map and its approximate.

The map f ρ (x) is class C 1 in x, and we assume that ∂fρ(x) ∂x

x is ultimately bounded for all x ∈ X . Thus, from Eq. ( 16), we write:

f ρ (x(k + j)) = f ρ (x(k + j -1)) + ξ ρ (k + j -1) (17)
+ ∂f ρ ∂x(k + j) x(k+j-1) f ∂ ρ (k+j-1)
δx(k + j -1) .

Expanding Eq. ( 17) along the horizon gives:

ρ(k + 1) = ρ(k) + f ∂ ρ (k)δx(k) + ξ ρ (k) , . . . ρ(k + N p -1) = ρ(k + N p -2) + ξ ρ (k + N p -2) + f ∂ ρ (k + N p -2)δx(k + N p -2) .
It is a fact that ρ(k) and δx(k) are known variables at each instant k, whereas f ∂ ρ (k) can be numerically evaluated. In practice,

f ∂ ρ (k +j) is unknown for j ∈ N [1,N -2]
. Then, in order to construct the estimation mechanism, we consider that

f ∂ ρ (k + j) = f ∂ ρ , ∀ j ∈ N [1,N -2]
, where f ∂ ρ denotes the partial derivative evaluated at instant k.

By doing so, and disregarding the residual terms, it is implied that ρ(k+j|k) = ρ(k+j-1|k-1)+f ∂ ρ δx(k+j-1|k). This is: the estimate for the future scheduling variable ρ(k +j|k) can be written as the sum of the estimate from the previous sample ρ(k + j -1|k -1) with a adjustment term f ∂ ρ δx(k + j -1|k). Accordingly, we can write the vector-wise extrapolation in a recursive fashion:

ρk [k+Np-1] = (ρ k-1 [k+Np-1] ) † + f ∂ ρ (δx [k+Np-2] ) † , ( 18
)
where the (•) † operator indicates a correction of the vectors with the known terms 10 . If sought, a forgetting factor can be added to (ρ k-1 [k+Np-1] ) † , multiplying the original vector by an exponentially decaying term, such as Ie -k/kmax . Such forgetting factor can be used to attenuate the amount of mistaken information passed from the scheduling sequence estimate ρk-1

[k+Np-1] to the following ρk [k+Np -1] . We note that (δx [k+Np-1] ) † can be computed thanks to the nominal prediction model of Eq. ( 9), using ρk-1 [k+Np-1] and v [k+Np-1] as the solution of the previous MPC iteration.

We stress that the recursive extrapolation law in Eq. ( 18) does not ensure that the scheduling sequence estimates abide to the scheduling parameter set P. Thus, each new extrapolation vector ρk

[k+Np-1] is "clipped" with respect to the bounds on ρ and δρ.

Lemma 2 Assume that input-to-state stability is ensured, despite the biased prediction model. Then, the recursive extrapolation algorithm in Eq. ( 18) converges. This is, after a finite amount of steps k c it holds that

lim k→kc ρk [k+Np-1] → ρ kc [k+Np-1] .
Proof 3 The convergence property can be demonstrated with the aid of the residual term ξ ρ (k + j), which should turn null. We demonstrate that this is true for j = 1; the same holds for the following steps due to equivalency. Considering that lim k→∞ x(k

+ 1) = x(k), take ξ ρ (k) = f ρ (x(k + 1)) -f ρ (x(k)) -f ∂ ρ δx(k).
Due to the stabilisation implied by the MPC, it directly follows that

lim k→∞ f ρ (x(k + 1)) = lim k→∞ f ρ (x(k)) and lim k→∞ δx(k) = 0. Then, lim k→∞ ξ ρ (k) = -lim k→∞ f ∂ ρ δx(k) → 0.
This concludes the proof.

Lemma 3

The estimation error (residual) is ultimately bounded to a convex set

E ⊂ P. It follows that ξ ρ ∞ = γ ρ + f ∂ ρ δx.
Proof 4 The residual term from Eq. ( 17) is:

ξ ρ (k + j + 1) = f ρ (x(k + j + 1)) -f ρ (x(k + j)) (19) -f ∂ ρ δx(k + j) .
Using a triangular inequality, we have:

ξ ρ (k + j + 1) ∞ ≤ f ρ (x(k + j + 1)) -f ρ (x(k + j)) ∞ + f ∂ ρ δx(k + j) ∞ . ( 20
)
As of Assumption 1, we have:

ξ ρ (k + j + 1) ∞ ≤ γ ρ δx(k + j) ∞ + f ∂ ρ δx(k + j) ∞ . Since f ∂ ρ is ultimately bounded, i.e. f ∂ ρ ∞ ≤ f ∂ ρ , it follows that: ξ ρ ∞ = γ ρ + f ∂ ρ δx. Since ξ ρ ∞ < ρ, E ⊂ P. This concludes the proof.
Remark 8 We note that even if f ρ (•) is not locally Lipschitz, different bounds for ξ ρ ∞ can still be computed, since ξ ρ converges to zero as of Lemma 2. The computation of these bounds in not the focus of this paper.

Robust NMPC for Tracking

In this Section, we detail the proposed robust NMPC scheme for Tracking. The proposed strategy holds similarities to the tracking NMPCs from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF][START_REF] Köhler | A nonlinear tracking model predictive control scheme for dynamic target signals[END_REF], but the method is conceived with regard to nominal predictions obtained through the qLPV model in Eq. [START_REF] Boccia | Stability and feasibility of state constrained MPC without stabilizing terminal constraints[END_REF]. The use of these nominal predictions makes our method faster than the prior, since the resulting numerical toughness is much closer to that of a QP. As done in [START_REF] Tito | A constraint-tightening approach to nonlinear model predictive control with chance constraints for stochastic systems[END_REF], we use constraint tightening, terminal cost and terminal constraints to ensure recursive feasibility and input-to-state stability properties.

Firstly, we make reference to Remark 9 from [START_REF] Köhler | A nonlinear tracking model predictive control scheme for dynamic target signals[END_REF]: the NMPC propositions in previous references [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF][START_REF] Köhler | A nonlinear tracking model predictive control scheme for dynamic target signals[END_REF] provides exponential closed-loop stability and recursive feasibility guarantees in the case of no model mismatch. Nevertheless, this is rarely the case in any practical application. When considering the qLPV predictions based on an estimated scheduling sequence, as we do in this paper, model mismatches inherently emerge. Therefore, in order to ensure robust recursive feasibility and stability (despite bounded disturbances and prediction mismatches), the MPC problem should be adjusted using constraint tightening techniques [START_REF] Tito | A constraint-tightening approach to nonlinear model predictive control with chance constraints for stochastic systems[END_REF]. This is what we perform in this paper with the use of the zonotopes detailed in Sec. 2.3. Through the sequel, we assume that robust constraint satisfaction is guaranteed thanks to the constraint tightening and the corresponding disturbance propagation zonotopes. Thence, we progress by detailing the proposed MPC setting for the case of possibly unreachable output reference signals.

In order to potentially increase the closed-loop domain of attraction and avoid loss of feasibility due to set-point changes, we benefit from the artificial referencing mechanism from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF]. Instead of tracking the time-varying output set-point y r we consider an additional decision variable y a . The tracking convergence is then ensured with the inclusion of an offset cost V O (y a -y r ), which penalizes the deviation between the artificial output reference y a and the real set-point y r .

The terminal set in ( 5) is adjusted to the tracking case, instead of the origin, considering the equilibria defined by any triad (x r , u r , y r ) which verifies [START_REF] Menezes Morato | Novel qLPV MPC design with least-squares scheduling prediction[END_REF]. Thus, we consider a terminal Tracking Robust Positive Invariant Set (TRPI set):

Definition 2 (Parameter-dependent TRPI Set) Consider a set Γ(ρ) ⊆ R nx+ny and a terminal control law u t = κ t (x, y r ). Γ(ρ) is a TRPI set for system (1) if, for all (x, y r ) ∈ Γ(ρ) and w ∈ S(N p ), (A(ρ)x + B(ρ)κ t (x, y r )) + w, y r ) ∈ Γ(ρ).

Remark 9

The TRPI sets differ from the Tracking Positive Invariant (TPI) sets from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF] due to the robustness properties. The synthesis of TPI sets implicitly considers that the real system trajectories and the nominal predictions are identical.

The previous definition implies that once a given state x and output reference y r are found inside the parameter-dependent set Γ(ρ), the terminal control law κ t : R nx+ny → R nu ensures that any subsequent state, with the same output reference, is also inside this set, regardless of the bounded load disturbance w ∈ W.

Thus, based on the current state measurement x(k), the estimated scheduling sequence ρk

[k+Np-1] , and the current output set-point y r , the NMPC optimisation ( 15) is re-stated as follows, being g x (•) and g u (•) the nonlinear maps that convert steady-state output references into state and input targets that satisfy Eq. ( 13), respectively:

min v,ya Np-1 j=0 (x(k + j|k) -x a , (v(k + j|k) + (21) K(ρ(k + j|k))x(k + j|k)) -u a ) +V (x(k + N p |k) -x a ) + V O (y a -y r ) , s.t. :                  x(k + j + 1|k) = A π (ρ(k + j|k))x(k + j|k) +B(ρ(k + j|k))v(k + j|k), j ∈ N [0,Np-1] , (x(k + j|k), v(k + j|k)) ∈ Z π (j), j ∈ N [0,Np-1] , x a = g x (y a ), u a = g u (y a ) (x(k + N p |k), y a ) ∈ Γ(ρ(k + N p -1)) .
In this problem, Z π (j) represent the tightened con-straints, as previously detailed, and Γ the terminal set. The set of states x(k) ∈ X a (N p ) such that this problem has a feasible solution is called the domain of attraction of the proposed controller 11 . We proceed by further detailing the MPC costs (•, •), V (•), and the terminal TRPI set Γ(•).

Assumption 3 (1) The stage cost function is positive definite and uniformly continuous such that:

(x, u) ≥ α ( x ) | (x 1 , u 1 ) -(x 2 , u 2 )| ≤ λ x ( x 1 -x 2 ) (22) + λ u ( u 1 -u 2 ),
where α , λ x and λ u are K-functions.

(2) The set of admissible artificial references

Y t = {y a ∈ R ny : (g x (y a ), y a ) ∈ Γ(g ρ (y a ))} is a convex subset of Y a (N p ) = {y a ∈ R ny : (g x (y a ), g v (y a )) ∈ Z π (N p )}. (3)
The offset cost function is positive definite, uniformly continuous and convex, thus assuring that the minimizer y o a = arg min ya∈Yt V O (y a -y r ) is unique, being V O (•) a weighted quadratic term. Furthermore, for any y r ∈ R ny and y a ∈ Y t , we have

V O (y a -y r ) -V O (y o a -y r ) ≥ α O ( y a -y o a )
, where α O is a K-function 12 .

Assumption 4 (1) The terminal control law satisfies κ t (x a , y a ) = u a for all admissible equilibrium points (x a , u a , y a ) 13 . (2) The terminal set Γ is an admissible TRPI set. That is, we have Γ as a subset of Λ(N p ) = {(x, y) ∈ R nx × Y t : (x, κ t (x, y)) ∈ Z(N p )} 14 satisfying Def. 2 for the control law u = κ t (x, y) and disturbances w ∈ S(N p ). (3) The terminal cost function V (x -x a ) is a Lyapunov function for the unconstrained qLPV system (1) such that for all (x, y a ) ∈ Γ(ρ) there exist constants b > 0 and σ > 1 such that V (x -x a ) ≤ b|xx a | σ . Thus, we have:

V (A(ρ)(x -x a ) + B(ρ)κ t (x - x a , y))-V (x-x a ) ≤ (x-x a , κ t (x, y)x-u a )
, where x a = g x (y a ) and u a = g u (y a ).

Notice that the previous assumptions presented are similar to the ones in [START_REF] Tito | A constraint-tightening approach to nonlinear model predictive control with chance constraints for stochastic systems[END_REF] and [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF], which develop control strategies for robust regulation and nominal tracking respectively, but herein extended to the robust case, with a 11 Notice that due to the freedom provided by the artificial reference ya, the feasibility property becomes independent of yr [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF]. 12 Notice that if the target output is admissible (yr ∈ Yt), the minimiser is y o a = yr and the previous inequality is reduced to VO(y) ≥ αO( y ). 13 An evident explicit alternative for this terminal law is κt(x, ya) = Kt(fρ(x))(x -gx(ya)) + gu(ya). 14 Note that Z(Np) is the set of points (x, u) such that u = v + K(fρ(x))x with (x, v) ∈ Zπ(Np).

qLPV prediction model. Accordingly, we provide a computationally elegant solution that can be used to compute parameter-dependent qLPV terminal ingredients, through the solution of some LMIs. The following Theorems provide recursive feasibility and exponential stability guarantees for the qLPV system in Eq. ( 1) subject to the MPC control law from Eq. ( 5). For the sake of presentation clarity, the proofs are provided in Appendix A.

Consider parameter-dependent TRPI set Γ(ρ) := X f (ρ) × Y t , being X f := x | x T P (ρ)x ≤ 1 an ellipsoidal terminal set. Furthermore, consider the sublevel terminal cost V (x) = x T P (ρ)x and the terminal qLPV feedback k t (x, y) = K t (ρ)(x -g x (y)) + g u (y). Note that the terminal scheduling parameter ρ is, by definition, ρ = f ρ (x). Since we seek x → g x (y), it is implicitly implied that the terminal scheduling parameter is ρ = g ρ (y). The NMPC optimisation can thus be implemented with the following terminal constraint: (x(k + N p |k), y a ) ∈ Γ(y a ) := Γ(f ρ (y a )), where y a is a free optimisation variable tuned through the offset cost V O to track the actual set-point y r .

In order to detach this nonlinear dependency, we synthesise ρ-dependent terminal ingredients, based on a positive-definite matrix P (ρ). Likewise, we use a set X f (ρ) which is robust positively invariant regarding the closed-loop dynamics (Eq (6)), in such way that Γ(ρ) is a TRPI set. Consider the tracking system e(k + j|k) = x(k + j|k)g x (y), with g x (y) (piece-wise) constant and x(k + j|k) being the nominal qLPV predictions from Eq. ( 9). We denote w † (k + j) := w(k + j) + ∆ j as the resulting (bounded) uncertainties (load disturbance summed to the model-process mismatch due to the differences between ρ(k + j|k) and ρ(k + j)). Considering the nominal law u(k) = k t (x, y), it follows that:

e(k + j + 1|k) = A π (ρ(k + j))e(k + j|k) + θ(k + j) ,(23) θ(k + j) = (A(ρ(k + j)) -I nx )g x (y) (24) + B(ρ(k + j))g u (y) + w † (k + j) ∈ Θ.
Accordingly, the following Theorems ensure that the error dynamics in Eq. ( 23) converge to the origin (i.e. tracking is ensured).

Theorem 2 from [9]

Consider Assumptions 3 and 4 hold. Suppose the there exists a terminal control law u = K t (ρ)(x -g x (y)) + g u (y). Consider that the MPC is given by Eq. (21), with a terminal state set given by X f (ρ) and a terminal cost V (x, ρ). Then, input-to-state stability is ensured if the following conditions hold ∀ρ ∈ P: (C1) The origin e = 0 lies in the interior of X f (ρ); (C2) Any consecutive error to e, in closed-loop given by e + = A π (ρ)e + θ lies within X f (ρ);

(C3) The discrete Lyapunov equation is verified within this invariant set, this is, ∀ x ∈ X f (ρ), ∀ ρ ∈ P, and ∀ δρ ∈ δP: V (e + , ρ + δρ) -V (e, ρ) ≤ -(x, k t (x, y)g u (y)). (C4) The image of the terminal feedback lies within the admissible control domain: K t (ρ)(x -g x (y)) + g u (y) ∈ U , ∀ρ ∈ P. (C5) The terminal set X f (ρ) is a subset of X .

Assuming that the initial solution of the MPC problem v [0,Np-1] is feasible, then, the MPC is recursively feasible, steering e = (x -g x (y)) to the origin.

Theorem 3 Tracking Robust Positive Invariant Set

Assume that there exists an ellipsoidal terminal set X f (ρ). X f is a robust positively invariant set iff, for any e ∈ X f and ρ ∈ P, i.e. e T P (ρ)e ≤ 1, it follows that (e + ) T P (ρ + δρ)e + ≤ 1, i.e. the successor state e + is also inside X f , which implies in:

(A π (ρ)e + θ) T P (ρ + δρ)(A π (ρ)e + θ) ≤ 1 . ( 25
)
Then, Γ(ρ) is a TRPI for system (1) as follows:

Γ := {(x, y) ∈ R nx×ny | (x -g x (y)) ∈ X f (ρ), (26) h(x, K t (ρ)(x -g x (y) + g u (y)) ∈ Y t } .
Theorem 4 Terminal Ingredients Conditions (C1)-(C5) of Theorem 2 and the inequality of Theorem 3 are satisfied if there exist a symmetric parameter-dependent positive definite matrix P (ρ) : R np → R nx×nx , a parameter-dependent rectangular matrix W (ρ) : R np → R nu×nx , and a scalar λ ∈ ]0, 1] such that Y (ρ) = (P (ρ)) -1 > 0, W (ρ) = K t (ρ)Y (ρ) such that LMIs (27)-( 28) and the BMI (29) hold for ρ + = ρ + δρ for all ρ ∈ P and δρ ∈ δP, and θ and w † given as the vertices of Θ and S(N p ), respectively, under the minimization of log det{Y (ρ)}.

Remark 10

The BMI in Theorem 4 can be solved through simple bisection search over the optimization plane since 0 < λ ≤ 1, by construction, as argues [START_REF] Yang | An optimal approach to output-feedback robust model predictive control of LPV systems with disturbances[END_REF].

Remark 11 Theorem 4 provides infinite-dimensional inequalities, which must hold ∀ ρ ∈ P and ∀ δρ ∈ δP. In practice, the solution can be found by enforcing the inequalities over a sufficiently dense grid of points (ρ, δρ) along the P × δP plane. Then, the solution can be verified over a denser grid. The parameter-dependency of P may be dropped if the system is quadratically stabilizable, but this may result in quite conservative performances.

Proposition 1 (Recursive Feasibility) Let there exist a solution Y (ρ) to Theorem 4. Then, given any x ∈ X a (N p ), y r ∈ R ny and u = κ(x, y r ), we have x + = A(f ρ (x))x + B(f ρ (x))u + w ∈ X a (N p ), ∀w ∈ W. Consider an optimal sequence v = (v 0 , v 1 , . . . , v Np-1 ) and an optimal artificial target y a (x, y r ) = y a . Then vc = (v 1 , . . . , v Np-1 , κ t (x(k +N p |k), y a )-u a ) and y c a = y a define feasible (candidate) solution of the MPC problem in Eq. (21) for any y r ∈ R ny and w ∈ W. This means that the NMPC optimisation in (21) is recursively feasible.

Proposition 2 (Error ISS) Let there exist a solution Y (ρ) to Theorem 4. Then, the qLPV system (1) in closed loop with the MPC input (5) has uniformly exponentially input-to-state stable error dynamics (as of Eq. ( 23)). That is, for any feasible initial condition x 0 and constant set-point y r ∈ R ny , with w(k) ∈ W, it is implied that:

x(k) -x a (k) ≤ β( x(0) , k) + γ(w), ( 30
)
where β and γ are respectively a KL-function and a Kfunction and w is such that w(k) ≤ w, ∀k.

Remark 12 Note that for an admissible equilibrium state x a the virtual control sequence v = (v a , . . . , v a ), where u a = v a + Kx a , with artificial reference y a is admissible, since it maintains the system at x a . Therefore, the set of admissible equilibrium states

X s (Z π (N p ), Y t ) = {x ∈ R nx : ∃ u a ∈ R nu , (x a , u a ) ∈ Z s (Z π (N p )), h(x a , u a ) ∈ Y t } is a subset of X a (N p )
and feasibility is not lost through set-point changes.

Remark 13

The feasible candidate (v c , y c a ) can be used as a starting point for (21) in order to reduce computational load.

Remark 14

In general, the functions g x (•) and g u (•) are not explicitly known. Therefore, rather than considering the constraints x a = g x (y a ) and u a = g u (y a ), we can include x a and u a as additional decision variables and enforce the analogous constraints x a = A(f ρ (x a ))x a + B(f ρ (x a ))u a , y a = C(f ρ (x a ))x a + D(f ρ (x a ))u a .

Artificial Reference Choice

The addition of the artificial reference y a ∈ Y t prevents feasibility loss due to changes in set-point or nonadmissible targets. However, it significantly increases the computational cost of the MPC optimization problem by introducing the new decision variable y a and associated nonlinear constraints x a = g x (y a ) and u a = g u (y a ).

In the sequel, a method for avoiding such additional computational burden is presented. First, consider the optimal admissible target y o a = arg min ya∈Yt V O (y a -y r ) and a feasible candidate artificial target15 y c a . Let y α a be a convex combination of y c a and y o a , i.e. y α a = (1 -α)y c a + αy o a , for an α ∈ [0, 1]. Note that, from the convexity of

                  Y (ρ) (A(ρ)Y (ρ) + B(ρ)W (ρ)) T Y (ρ) W T (ρ) 0 (A(ρ)Y (ρ) + B(ρ)W (ρ)) T 0 0 (A(ρ)Y (ρ) + B(ρ)W (ρ)) Y (ρ + ) 0 0 θ Y (ρ + ) 0 0 Y (ρ) 0 Q -1 0 0 0 In x 0 W (ρ) 0 0 R -1 0 0 0 In x 0 θ T 0 0 0 θ T 0 0 (A(ρ)Y (ρ) + B(ρ)W (ρ)) Y (ρ + ) 0 0 θ Y (ρ + ) 0 0 0 0 In x 0 0 0 In x 0 0 0 0 In x 0 0 0 In x                   ≥ 0 , (27) (a)
  ui -I nu,{i} ur 2 I nu,{i} W (ρ) (I nu,{i} W (ρ)) T Y (ρ)   ≥ 0 , ∀i ∈ N [1,nu] , (b)   xi -I nx,{j} xr 2 I nx,{j} Y (ρ) (I nx,{j} Y (ρ)) T Y (ρ)   ≥ 0 , ∀j ∈ N [1,nx] , (28)
       λY (ρ) (A(ρ)Y (ρ) + B(ρ)W (ρ)) T 0 (A(ρ)Y (ρ) + B(ρ)W (ρ)) T (A(ρ)Y (ρ) + B(ρ)W (ρ)) Y (ρ + ) θ Y (ρ + ) 0 θ T (1 -λ) θ T (A(ρ)Y (ρ) + B(ρ)W (ρ)) Y (ρ + ) θ Y (ρ + )        > 0 . ( 29
)
Y t , y α a ∈ Y t . Thus, consider the following auxiliary optimization problem:

max v,α α (31) s.t.
:

                 x(k + j + 1|k) = A π (ρ(k + j|k))x(k + j|k) +B(ρ(k + j|k))v(k + j|k), j ∈ N [0,Np-1] , (x(k + j|k), v(k + j|k)) ∈ Z π (j), j ∈ N [0,Np-1] , α ∈ [0, 1] (x k+Np|k , y α a ) ∈ Γ(y α a )
, Based on the solution α * of (31), we can then solve the simpler optimization problem [START_REF] Tito | A constraint-tightening approach to nonlinear model predictive control with chance constraints for stochastic systems[END_REF], with y r = (1α * )y c a + α * y o a , to obtain the optimal control sequence v [k,k+Np-1] . This two-step optimization is analogous to (21) with an arbitrarily large offset cost V O (•) (prioritizing the proximity of the artificial target from the optimum y o a ) and the added restriction that y a must be inside the line segment [y c a , y o a ].

Remark 15 Note that if, for some time instant k, the solution of (31) is α * = 1, we will have, assuming a constant reference, y c a = y o a for the future instants. The auxiliary problem (31) thus becomes irrelevant until there is a change on y r .

Remark 16 Since α is a real scalar in [0, 1], an alternative to the optimisation in Eq. (31) is to perform a bi-section 16 mechanism with α c k ∈ [0, 1]. By doing so, we test the feasibility of (15) for a given α c k . If the problem is feasible, we apply the correspondent control law; otherwise, we take α = 0 (which is by definition feasible). A viable choice of candidates α c k , in order to ensure the convergence of y a to y o a , can be of

α c k = 1 if α c k-1 was feasible and α c k = 1 2 α c k-1 if not (bisection)
. This reduces the computational complexity to a maximum of two optimization problems given by [START_REF] Tito | A constraint-tightening approach to nonlinear model predictive control with chance constraints for stochastic systems[END_REF], without the undesired nonlinearities.

Results

Application and qLPV Embedding

In this Section, we use a numeric 17 benchmark from the literature to demonstrate the properties and effectiveness of our method. We adapt the cascaded tank process from [START_REF] Henrik | The Quadruple-tank Process: A Multivariable Laboratory Process with an Adjustable Zero[END_REF], considering the interconnection of two tanks, with an open hole at the bottom of the first tank, which leaks fluid to the second. The latter has a pump at its end, regulated by a local proportional controller. The resulting nonlinear level dynamics are:

dh 1 (t) dt = - a 2gh 1 (t) A 1 + γ A 1 u(t) , dh 2 (t) dt = a 2gh 1 (t) A 2 - k A 2 h 2 (t) .
Each h i (t) represents the water level at the i-th tank; u(t) represents the tension applied for the main pump, for which the corresponding flow is γu(t). We are mainly interested with the dynamics of y(t) = h 1 (t). The tank cross sections A i are of 1 cm 2 , while the outlet hole cross section a is of 0.05 cm 2 . The pump parameter γ is of 1.4 cm 3 /Vs. The proportional coefficient k is of 1.1. The process constraints are:

h j ∈ [1, 10] cm, ∀j ∈ N [1,2] (level, process variable); u ∈ [0, 5] V (pump flow, control input).
This nonlinear model is Euler-discretized with T s = 250 ms, and qLPV-embedded to Eq. ( 1) with:

A(ρ) = (1 -Tsa √ 2g A1 ρ) 0 Tsa √ 2g A2 ρ (1 -Tsk A2 ) , B(ρ) = Tsγ A1 0 , ρ = f ρ (x) = (x 1 ) -0.5 ∈ [0.3, 1] cm -0.5 .
In order to compute the bounds on δρ, one can use either constrained optimisation or interval arithmetic methods as described in Remark 3. The obtained bounds are defined through δρ := δρ ∈ R 1 : -0.034 ≤ δρ ≤ 0.0052 . The process is subject to bounded additive uncertainties w(k) ∞ ≤ 0.05 cm. In the sequel, we compare the proposed robust qLPV NMPC with the NMPC for Tracking algorithm from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF].

Terminal Ingredients

Consider the convex output tracking set Y t := [1 , 10] cm. The output steady-state condition from Eq. ( 13)

implicitly defines g x (y r ) = y r a √ 2gyr k T and g u (y r ) = a √ 2gyr γ
. We use a prediction horizon of N p = 4 steps and the quadratic stage cost (x, u)

= x 2 Q + u 2 R with Q = I nx and R = 1.
In order to synthesise the terminal ingredients via Theorem 4, we partition the output set Y t into ten partitions, thus finding one parameter-dependent RPI set Γ(ρ, y r ), for which y r parametrises the partitions, i.e. y r ∈ [1, 2] or [START_REF] Victor | Nonlinear robust predictive control with distrubance propagation via zonotopes[END_REF][START_REF] Cisneros | Nonlinear model predictive control for models in quasi-linear parameter varying form[END_REF] or [START_REF] Cisneros | Nonlinear model predictive control for models in quasi-linear parameter varying form[END_REF][START_REF] Hoffmann | A survey of linear parameter-varying control applications validated by experiments or high-fidelity simulations[END_REF], and so on up to [START_REF] David Q Mayne | Constrained model predictive control: Stability and optimality[END_REF][START_REF] Mohammadpour | Control of linear parameter varying systems with applications[END_REF]. For the NMPC algorithm, we use terminal ingredients synthesised through the procedure from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF]Appendix B].

Figure 1 shows the parameter-dependent TRPI sets (for frozen values of ρ ∈ P) used for the proposed algorithm and the quadratic TPI sets used for the NMPC algorithm [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF]. All sets are translated to from the error coordinates (x -x r ) to the state coordinates x, centered at the different state targets x r of each partition. As we can see, Theo. 4 generates sufficiently a large terminal region, meaning that any reference point y r ∈ Y t can be tracked.

Scheduling Sequence Estimation

Before presenting the actual control results, we provide the scheduling sequence extrapolation estimates obtained with the recursive method presented in Sec. 3. As detailed in Lemma 2, convergence is indeed verified. The obtained bounds for the estimation error are ξ ρ ∞ ≤ 0.015 cm -0.5 , as shown in Fig. 2. The extrapolation mechanism offers very precise estimates ρk [k+Np -1] , which means that the nominal qLPV predictions obtained through Eq. ( 9) are very close to the real system trajectories of Eq. ( 8) and thus the disturbance propagation along the horizon is reduced.

Disturbance Propagation

Now, we show the disturbance propagation reachable sets S(j), ∀j ∈ N [0,Np] for the proposed algorithm, considering the zonotope disturbance propagation paradigm. Note that S(0) stands for the load disturbance set W; the following zonotopes comprise the propagation of the load disturbances and the modelprocess mismatches along the horizon. Thse sets are computed according to Theorem 1. In Fig. 3, we show the collection of sets S(j) over the x 1 × x 2 plane (Def. 1). We recall that the vertices of S(N p ) were used to construct the terminal ingredients through Theorem 4.

Simulation Scenarios

We consider two different simulation scenarios. Since the NMPC algorithm from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF] is not robust by design, we first consider and compare the obtained tracking performances of both algorithms, without the presence of load disturbances (i.e. with w nil). Then, we consider the disturbance rejection robust performances solely of the proposed robust algorithm.

Nominal Tracking Performances

Considering a step-like piece-wise constant output target signal which passes through y r = 9, 5, 4 and 2 cm, the obtained tracking performances with both algorithms are shown in Fig. 4: sub-figure (a) presents the resulting state, input and output trajectories, while sub-figure (b) shows the state phase plane and the terminal sets. Complementary, Fig. 5 provides the values for the artificial reference tuning variable α for the proposed mechanism.

As we can see, the obtained tracking performances both methods are offset-free steady-state output points. The proposed method ensures slightly faster convergence than the original NMPC for Tracking scheme from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF].

The main advantage resides in its simpler implementation, of QP-alike numerical burden, enable through the qLPV embedding. With the qLPV model, nonlinear mapping do not have to be solved internally by the optimisation procedure. The proposed qLPV NMPC mechanism requires only the operation of: one linear recursive law (Eq. ( 18)) and one QP problem (Eq. ( 21)). In the moments of reference changes, one bisection search (Remark 16) is also required, which increases the number of QPs to at most 5 iterations. In contrast, the original NMPC for Tracking requires the solution of an NPhard optimisation problem per sampling instant, which is numerical-wise much harder.

In order to better compare the two tracking controllers, we assess the obtained performance results with performance indexes, presented in Table 1. We discuss them:

• Firstly, we note that there is an overall performance enhancement: there is a small decrease on the integral of absolute output tracking error (IAE) w.r.t.

[8], of roughly 12 %. • Complementary, we stress that the generated control input is smoother with the proposed method: we obtain 43 % reduced total variance (TV) of the control signal. • Furthermore, the average computational time needed to solve the control problem (t c ) is over 44 % smaller with the proposed qLPV method. This is very significant, since this is a small system (n x = 2). The complexity of the NP solution from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF] grows exponentially with (N p × n x ), which can a serious issue with time-critical systems. The proposed method has QP-alike burden, thus t c grows linearly with (N p × n x ), which means it is readily-conceived for embedded applications. • Finally, note that the model-process discrepancies (differences between Eqs. ( 8) and ( 9)) are very wellhandled with the zonotopic constraint-tightening approach, since the generated sets S(j) are small (Fig. 3). This corroborates previous discussions in the literature indicating this option as a promising alternative, see e.g. [START_REF] Köhler | A nonlinear tracking model predictive control scheme for dynamic target signals[END_REF]. We also stress that the terminal ingredients conceived with the proposed Theorem generate sufficiently large terminal sets, able to guarantee recursive feasibility for a rather larger output set.

Robust Tracking Performances

In order to illustrate the robustness properties of the proposed algorithm, we provide a second simulation scenario 18 . Consider three uniformly random disturbance sequences with unitary seeds, multiplied by decaying exponential terms, as illustrated in Fig. 6.

In Fig. 7, we show the state and output behaviours w.r.t. a step-like output target goal y r . Clearly, robust stability is ensured: as the load disturbance sequences dissipate, the error trajectories (x -x r ) converge to the origin; moreover, while w(k) is non-null, the states stabilize at constant steady-states regimes, as close as possible to x r . The average tracking error when disturbances are present is of 5.68 %, which is arguably small.

Conclusion

In this paper, we provided a novel Tracking NMPC algorithm through qLPV embedding. The method is based on the recursive estimation of the future qLPV scheduling trajectories, made available through a simple Taylor expansion. The propagation of the model-process prediction mismatches along the NMPC horizon are bounded with zonotopes. Furthermore, we provide an LMI-solvable remedy for the case of bounded additive disturbances, which computes a robust LPV feedback gain and parameter-dependent terminal ingredients. The derived tracking robust positive invariant set ensures recursive feasibility of the optimisation procedure as well as input-to-state stability of the process. We compare our method against the nominal tracking NMPC framework from [START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF], considering a cascaded tanks benchmark system. We find very similar tracking performances with much smaller computational stress with the proposed scheme, which benefits from the linear predictions of the qLPV realisation. The method is ready for embedded applications (the online stress is similar to that of a QP) and offers robustness towards bounded load disturbances with reduced conservatism.

with M ≥ 0 and:

M 11 = Y (ρ) -Y T (ρ)QY (ρ) -W T (ρ)RW (ρ) (A.7) -Y T (ρ)A T π (ρ)P (ρ + )A π (ρ)Y (ρ) , M 12 = -Y (ρ)A T π (ρ)P (ρ + )θ , (A.8) M 21 = -θ T P (ρ + )A π (ρ)Y (ρ) ,
(A.9) M 22 = -θ T P (ρ + )θ .

(A.10)

Using e = eP (ρ)Y (ρ) and W (ρ) = K t (ρ)Y (ρ) leads to:

(A π (ρ)e + θ) T P (ρ + )(A π (ρ)e + θ) -e T P (ρ)e (A.11) +e T Qe + e T K T (ρ)RK t (ρ)e ≤ 0 .
This inequality is a sufficient condition for (C3) with V (•) as a sub-level of X f . The fourth and fifth conditions (C4-C5) are verified by the direct application of the Schur complement to Eq. (28a) and Eq. (28b), respectively, using W (ρ) = K t (ρ)Y (ρ). They lead, respectively, to: u i -I nu,{i} u r 2 ≥ (A.12)

I nu,{i} K t (ρ) Y (ρ)P (ρ)Y (ρ) I nu,{i} K t (ρ) T .

Since the maximum normed F e of an e that belongs to some ellipsoid e T P e ≤ 1 is given by F (P -1 ) F T , it holds that the first inequality implies that the projection I nu,{i} K t (ρ)e (i.e. i-th control signal) is upper-bounded, in norm, by u i -I nu,{i} u r , which satisfies (C4). Analogously, LMI (28b) ensures that the projection I nx,{j} x (i.e. j-th state) is norm-bounded by x j -I nx,{j} x r , which satisfies condition (C5). This concludes the proof.

A.4 Proof of Proposition 1 (Recursive Feasibility)

Let Assumptions 3-4 hold. Consider there exists a solution Y (ρ) to Theorem 4. Then, from conditions (C1), (C2), (C4), and (C5) from Theorem 2, we can infer that the generated control signal provides recursively feasible solutions to the NMPC optimisation in (21). Take κ(x, y r ) = K t (ρ)(x -g x (y r )) + g u (y r ) with ρ = f ρ (x) (implied implicitly through the qLPV prediction model).

Then, the if the initial condition x(0) generates a feasible v, all future iterations of the optimisation will also be feasible: the generated control law control is admissible (condition (C4)) and all state variables x ∈ X f generate successor state variables x + which are also inside X f (condition (C2)), which contains the origin (condition (C1), terminal condition for (x -x r ) = (x -g x (y r )) and is a sub-set of X (condition (C5), which ensures that x, x + are admissible). This concludes the proof.

A.5 Proof of Propostion 2 (Error ISS)

Let there be a terminal stage cost V (•) such that Assumption 3 holds. Let Assumption 4 also hold and Theorem 1 be verified. Note that since (x -x a , u -u a ) is a quadratic stage cost x -x a 2 Q + u -u a 2 R , α , γ x and γ u indeed exists. Consider there exists a solution Y (ρ) to Theorem 4. Then, the closed-loop is stable due to (C3) of Theorem 2, which conversely ensures that δV (k) = V (x(k) -x a (k))

-V (x(k -1) -x a (k -1)) ≤ -x(k) -x a (k) Q + γ V (w) .

Assume that lim k→+∞ y a (k) → y o a . Analogously, use lim k→+∞ x a (k) → x o a := g x (y o a ). Then, thanks to the error dynamics in Eq. ( 23), we obtain:

x(k) -x a (k) Q ≤ β( (x(0) -x a (0)) , k) + γ(w) .

Since Q > 0 (and positive definite) and x a (0) = 0 (by design), we have x(k) -x a (k) Q ≥ x(k) -x a (k) and thus error-ISS is established. Now we consider the convergence of V O (y a (k) -y r ) such that the limit lim k→+∞ y a (k) → y o a holds. Let us define ŷa = (1 -α)y a (k) + αy o a , where α ∈ [0, 1] is the optimal solution from Eq. (31). From the convexity of V O (•), we obtain:

V O (ŷ a -y r ) ≤ (1 -α)V O (y a (k) -y r ) + αV O (y o a -y r ) .

We can use the Lipschitz continuity of the map x r := g x (y r ) in order to obtain x a (k)-xa ≤ L x y a (k)-ŷa , where L x > 0 is the Lipschitz constant of g x (•). Consider (y a (k) -ŷa ) = α(y a (k) -y o a ).

Since the closed-loop is stable, it follows that the total MPC cost dissipates over time, which implies in: V O (y a (k) -y r ) → γ y (w), we obtain V O (y a (k) -y r ) -V O (y o a -y r ) ≤ a 1 L σ x α σ-1 y a (k) -y o a σ + γ n (w), with σ > 1 and a 1 > 0 as a constant scalars. Taking the limit at both sides of this inequality leads to: .

Note that, in nominal conditions (reachable reference y r ∈ Y and null disturbances), we obtain γ Vo (w) = 0 and y o a = y r , and thus lim k→+∞ V O (y a (k) -y r ) → 0, which means the steady-state target is reached. In the case the reference isn't reachable and there are disturbances, we can only infer that lim k→+∞ y a (k) exists within (-µ(y o a , y r , w), +µ(y o a , y r , w)). Nevertheless, it is implied that lim k→+∞ V O (y a (k) -y r ) is bounded, which means it converges and error-ISS holds. This concludes the proof.

Fig. 1 .

 1 Fig. 1. Synthesised TRPI set partitions Γ(ρ) and the quadratic TPI set partitions.

Fig. 2 .

 2 Fig. 2. Scheduling Sequence Estimates ρk [k+Np-1] at different samples k and estimation error ξρ (and bounds).

Fig. 3 .

 3 Fig. 3. Zonotopic sets S(j).

  -Limón et al., 2018 Trajectories Quadratic TPI Set Parameter-dependent TRPI Set Proposed qLPV NMPC for Tracking Trajectories (b) State Phase Plane and TRPI sets.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Nominal Performances: Proposed qLPV NMPC vs. NMPC for Tracking [8].

Fig. 7 .

 7 Fig. 6. Disturbance Scenario.

V

  O (y a (k) -y r ) -V O (ŷ a -y r ) ≤ V (x a (k) -xa ) + γ Vo (w) ≤ a 1 x a (k) -xa σ + γ Vo (w) ≤ a 1 (L x y a (k) -ŷa ) σ +γ Vo (w) ≤ a 1 L σ x α σ y a (k) -y o a σ + γ Vo (w) .Then, from the convexity of V O (•), we have:V O (ŷ a -y r ) ≤ (1 -α)V O (y a (k) -y r ) + αV O (y o a -y r ) ,Thus, since 1 α > 1, we obtain:V O (y a (k) -y r ) ≤ a 1 L σ x α σ-1 y a (k) -y o a σ + γ Vo (w) + 1 -α α V O (y a (k) -y r ) + V O (y o a -y r ) .

Finally, since lim k→+∞ 1

 1

  V O lim k→+∞ y a (k) -y r ≤ V O (y o a -y r ) + γ n (w) . Note that V O (•) is a weighted quadratic cost by definition, thus | y a -y r | ≤ y a -y r T ≤ V O (y a -y r ). Thus: lim k→+∞ y a (k) -y r ≤ y o a -y r T + γ n (w) , -( y o a -y r T + γ n (w)) + y r ≤ lim k→+∞ y a (k) ≤ y o a -y r T + γ n (w) + y r µ(y o a ,yr,w)

Table 1

 1

	Performance Comparison.
						Method		IAE := k yr(k) -y(k)	TV := k u(k + 1) -u(k)	tc
						NMPC for Tracking [8]			108.23	35.93	90.53 ms
					Proposed qLPV NMPC for Tracking		95.31	20.44	50.55 ms
	10
								y r
	1 (cm)	6 8				NMPC for Tracking -Limón et al., 2018 [13] Proposed qLPV NMPC for Tracking
	x
			4
			2
				0	50	100	150	200	250	300
	2 (cm)	4 6
	x
		2
	u (V)	0 5	0 0	50 50	100 100	150 150	200 200	250 250	300 300
							Time (s)
					(a) State, Output and Input Trajectories.

We stress that the application results provided in[START_REF] Köhler | A nonlinear tracking model predictive control scheme for dynamic target signals[END_REF] solve the online NPs very fast, but relying on a solver-based solution (CaSaDi), which approximates the solution of the NPs.

[START_REF] Victor | Nonlinear robust predictive control with distrubance propagation via zonotopes[END_REF] The proposed method also considers bounded process disturbances. The uncertainty propagation zonotopes offer a direct extension for such case.

The bounds xj and uj can be found from Eq. (2) through linear programming.

It follows that P ⊆ fρ(X ), as detailed in[START_REF] Cisneros | Nonlinear model predictive control for models in quasi-linear parameter varying form[END_REF].

For such, the system must be ultimately bounded.

In this case, K must be parameter-independent such that Aπ(ρ) = A(ρ) + B(ρ)K becomes affine.

Bounds can be found from Eqs. (1)-(2) through interval arithmetic.

Note that at, instant k, the scheduling variable ρ(k) and the state deviation δx(k) are known. Thus, (ρ k-1[k+Np-1] ) † stands for ρk-1[k+Np-1] with the term ρ(k|k -1) replaced by ρ(k), whereas (δx [k+Np-1] ) † stands for δx [k+Np-1] with the term δx(k|k -1) replaced by δx(k).

From the Recursive Feasibility property (Theorem 1), the previous artificial reference y * a (x k-1 , yr) can be used as the candidate y c a .

In order to preserve the ISS property, this bisection search should also ensure that the terminal offset cost decays with the new artificial target, this is: V (x(k + 1|k) -gx(y α a )) -V (x(k|k) -gx(y c a )) ≤ -x(k|k) -gx(y c a ) Q.

2.4 GHz, 8 GB RAM Macintosh computer.

We opt not to test the method from[START_REF] Limon | Nonlinear MPC for tracking piece-wise constant reference signals[END_REF] against load disturbances since it is not a robust algorithm, which would in turn result in an unfair comparison.

Acknowledgements

This work has been supported by CNPq (304032/2019-0), CAPES (001) and ITEA3 European project (15016) EMPHYSIS. The Authors declare no conflict of interest.

A Proofs of Theorems and Propostions

A.1 Proof of Theorem 2

This proof is standard and basically gathers the sufficient conditions from Assumptions 3 and 4. Refer to [START_REF] David Q Mayne | Constrained model predictive control: Stability and optimality[END_REF].

A.2 Proof of Theorem 3

It is trivial to verify Eq. (25) by using x(k), x(k + 1) ∈ X f , ∀k ∈ N and y = h(x, K t (ρ)(x-g x (y))+g u (y)) ∈ Y t .

A.3 Proof of Theorem 4

We begin by showing the positive invariance of the ellipsoid. Applying the S-procedure, with λ > 0 to (25) and 1 -e T P (ρ)e ≥ 1 , we get:

which can be rewritten as:

with N > 0 and:

Applying a Schur complement over P (ρ+δ) for each N ij leads to (29). This ensures Theorem 3.

Complementary, we proceed by demonstrating that the resulting P (ρ) satisfies all five conditions of Theorem 2. (C1) trivially holds due to the ellipsoidal form of X f . (C2) is verified due to the fact that X f is a sub-level set of the terminal cost V (•). Therefore, if condition (C3) is verified, (C2) is consequently ensured.

The discrete Ricatti condition (C3) is verified through the solution of LMI (27). Since Q -1 > 0, R -1 > 0 and Y (ρ + δρ) > 0, we apply two consecutive Schur, complements, which leads to:


Aurélien Duret
email: aurelien.duret@ifsttar.fr

Yufei Yuan
email: yuan@tudelft.nl

Yuan Duret

Traffic State Estimation Based on Eulerian and Lagrangian Observations in a Mesoscopic Modeling Framework

come

INTRODUCTION

State of the art

Traffic state estimation (TSE) is crucial in real-time dynamic traffic management and information applications. The essence of TSE is to reproduce traffic conditions based on available observation data. One class of available estimation methods does not make use of traffic flow dynamics, but relies on basic statistics and interpolation. These are referred to as data-driven methods. Another class of estimation methods relies on dynamic traffic flow models. These are referred to as modelbased methods. The focus of this article is on the latter because it potentially provides better results than the former class in non-recurrent situations (work zones, accidents, social events, etc.), regarding the monitoring-forecasting capabilities.

Model-based TSE relies on two components: a model-based component and a data assimilation algorithm. The model-based component consists of two parts: (i) a dynamic traffic flow model to predict the evolution of the state variables; and (ii) a set of observation equations relating sensor observations to the system state. Thereafter, a data-assimilation technique is adopted to combine the model predictions with the sensor observations. For example, the Kalman filter (KF) [START_REF] Herrera | Incorporation of lagrangian measurements in freeway traffic state estimation[END_REF] and it advanced relatives, such as Extended KF [START_REF] Wang | Real-time freeway traffic state estimation based on extended kalman filter: A general approach[END_REF], Unscented KF [START_REF] Ngoduy | Applicable filtering framework for online multiclass freeway network estimation[END_REF], Ensemble KF [START_REF] Work | An ensemble kalman filtering approach to highway traffic estimation using gps enabled mobile devices[END_REF] have been widely applied in the field of traffic state estimation.

The same traffic flow model can be formulated in three two-dimensional coordinates regarding space x, time t and vehicle number n. Laval and Leclercq [START_REF] Laval | The hamilton-jacobi partial differential equation and the three representations of traffic flow[END_REF] have presented three equivalent variational formulations of the first-order traffic flow models, namely N (x, t) model, X(t, n) model, T (n, x) model respectively, under the theory of Hamilton-Jacobi partial differential equations. Under such defined coordinate systems, sensor observations from road networks can be defined into two categories: (i) Eulerian sensing data -observations (e.g., aggregated speeds, flows) from spatially-fixed sensors (such as inductive loops, video sensors, and radar sensors) over a fixed report frequency, this type is dominating the information sources in the field of transportation research for decades; and (ii) Lagrangian sensing data -information from probe samples at a fixed time interval (such as position and speed information of individual vehicles [START_REF] Herrera | Incorporation of lagrangian measurements in freeway traffic state estimation[END_REF], and/or probe spacing and position information [START_REF] Seo | Probe vehicle-based traffic state estimation method with spacing information and conservation law[END_REF]), this class is becoming an increasingly popular source. In literature, most of TSE applications are based on the traditional space-time (Eulerian) formulation. Aggregated traffic quantities (e.g., flows, densities or speeds [START_REF] Herrera | Incorporation of lagrangian measurements in freeway traffic state estimation[END_REF][START_REF] Wang | Real-time freeway traffic state estimation based on extended kalman filter: A general approach[END_REF][START_REF] Work | An ensemble kalman filtering approach to highway traffic estimation using gps enabled mobile devices[END_REF]) are usually considered as system states, but no individual vehicle tracking is involved. The popularity of this formulation is due to the fact that incorporating Eulerian data is straightforward and intuitive. Recent studies have shown that a first-order (LWR) traffic flow model can be formulated and solved more efficiently and accurately in vehicle number-time (Lagrangian-time) coordinates [START_REF] Leclercq | The lagrangian coordinates and what it means for first order traffic flow models[END_REF]. And its related Lagrangian formulation of TSE enables more accurate and efficient application of data assimilation methods, due to the solution to the mode-switching problem (traffic information travels in one direction), less non-linearity of the system model, and the nature set of observation equations to deal with Lagrangian data [START_REF] Yuan | Real-time lagrangian traffic state estimator for freeways[END_REF]. However, the computation cost depends on the discretized platoon size (set to 1 vehicle classically) and time grid (often set around 1 second), which might be time consuming. To the best of our knowledge, none of previous research has provided a complete TSE framework for assimilating both Eulerian and Lagrangian data under a vehicle number -space (Lagrangian-space:= L-S) formulation.

Objectives and contributions

This paper presents a generic data assimilation framework based on a mesoscopic-LWR model formulated in Lagrangian-space coordinates, using both Lagrangian and Eulerian observations. The term mesoscopic is in response to the two other counterparts, since the Lagrangian-time coordinates can apply in a microscopic simulation framework and the Eulerian coordinates can accommodate in a macroscopic one. In this work, the system model is the Lagrangian-space formulation of the LWR model. It individually represents vehicles but only tracks their states at cell boundaries. We will develop algorithms and observation models to incorporate data from both Eulerian and Lagrangian sensors. And we do not apply specific data assimilation techniques; instead we try to demonstrate the sequential data assimilation concepts via reasonable assumptions. The algorithms on how to estimate network traffic states under the proposed model-based framework from the two data sources will be the main contribution of this work.

Contents of the paper

This paper is organized as follows. Section 3 presents the underlying traffic dynamics, including its formulation, solutions and properties. Section 4 describes the methodology of the proposed TSE framework, including how to assimilate Lagrangian data and combine with Eulerian data. Sections 5 and 6 illustrate the model validation and an application to a freeway corridor. Discussion and conclusions are drawn in Sections 7 and 8.

LS-LWR MODEL

This section defines the underlying process model in the state estimation framework, where the model formulation, numerical solution and its properties are discussed.

Conservation law and variational theory

This section first presents a mesoscopic formulation of the LWR model as the process model in the estimation framework. The LWR model is formulated in vehicle platoon and space (n, x) coordinates. The current mesoscopic formulation combines a vehicular description with macroscopic behavioral rules. It relaxes the temporal coordinate, and this entitles a transformation of a temporal progressing approach (e.g., in Eulerian or Lagrangian-time simulation framework) to an event-progressing approach (trigger event can be the change of time headway or pace, and/or a correction procedure based on an observation from fixed loops or probe vehicles).

The formulation follows the principle of the Hamilton-Jacobi (HJ) theory, to find an expression of the LWR model in Lagrangian-space coordinates. This model is also referred to as the T -model.

The LWR model can then be described by a hyperbolic equation:

∂ x h -∂ N (1/V (h)) = 0 (1)
Here, h denotes the time headway. The inverse speed 1/v (or called pace τ ) can be derived from the fundamental diagram 1/V (h).

Previous authors have proposed to apply variational theory in Eulerian coordinates (x, t)

(1) and Lagrangian coordinates (n, t) [START_REF] Leclercq | The lagrangian coordinates and what it means for first order traffic flow models[END_REF]. Here, we transpose the demonstration in Lagrangianspace coordinates (n, x), following the same rationale in [START_REF] Leclercq | The lagrangian coordinates and what it means for first order traffic flow models[END_REF]. The problem can be expressed as the Hamilton-Jacobi derived from the fundamental diagram:

∂ x T = 1 V (∂ N T ) (2)
Here, the function 1/V represents the flux function of the problem.

Numerical solutions in Lagrangian-space coordinates

Here, a Godunov scheme is applied to solve the conservation law equation (hyperbolic equation)

above with an upwind method. This would preserve the numerical benefit of Lagrangian traffic flow models. Fig. 1 illustrates the mesoscopic numerical grid (see grey area) for the Godunov scheme. On the mesoscopic grid, the time headway is determined by:

h x+∆x i = h x i + ∆x ∆n .( 1
V (h x i ) - 1 V (h x i-1 ) ) (3)
The CFL condition that guarantees the convergence of the Godunov scheme is:

∆n ≥ max h ∂ h ( 1 V ) ∆x (4)
Alternatively, the problem can also be expressed in terms of T (n, x) considering the 'passage time' flux that crosses the boundary of the cell n, regarded as a variational formulation of the T -model:

T (n, x) -T (n, 0) ∆x = 1 V ( T (n,x)-T (0,x) ∆n ) (5)
In this expression, V depends on the fundamental diagram. Here, we consider a triangular fundamental with three parameters: the free-flow speed v m , the maximum wave speed w and the jam density k x . It can be expressed by:

1 V (h) = max( 1 v m , -k x ( 1 wk x -h)) (6)
The numerical solution to the problem is simplified as:

T (n, x) = max(T (n, 0) + x v m , T (0, x + n k x ) + n w.k x ) (7)
Finally, the origin n = 0, x = 0 could be shifted to n -∆n, x -∆x and we find:

T (n, x) = max(T d , T s ), (8)
where T d = T (n, x -∆x) + ∆x vm and T s = T (n -∆n, x -∆n kx ) + ∆n w.kx represents the demand and the supply term respectively. The demand term defines the arrival time of a vehicle from upstream in non-constrained (free-flowing) conditions. The resulting passage time of a vehicle is at least equal to its arrival time but could be delayed due to the downstream conditions. Thus, the supply time provides such information in constrained (congested) conditions.

This numerical solution indicates traffic flow is divided into vehicle platoons of certain size ∆n, and road stretch is discretized into spatial cells of certain length ∆x. Note that, the cell length ∆x in the simulation is not necessarily to be equal. The state in this formulation is the passage time T (n, x) of vehicle platoons at cell boundaries. This state is always determined by the maximum of two uncorrelated terms: the demand (arrival) time and the supply time. For an elaborate description we refer to [START_REF] Laval | The hamilton-jacobi partial differential equation and the three representations of traffic flow[END_REF].

Properties

The current mesoscopic formulation is based on the notions from the variational theory. It can incorporate the numerical benefits and modeling flexibility of both Eulerian and Lagrangian-time models. Simultaneously, this formulation allows state distinction on both cell class and vehicle class, combining a vehicular description with macroscopic behavioral rules. It individually represents vehicles (platoons) but only tracks their passage times at cell boundaries. Therefore, travel times can be easily derived from the model, which is more convenient compared to other (e.g., Eulerian or Lagrangian) formulations of state estimation. This discrete model evolves state by state, with only one expression to consider all traffic conditions. Hence, it does not require memory and it is more flexible and time-efficient for data assimilation (no complex matrix inversion and multiplication). Moreover, the numerical scheme allows for long cells and cell boundaries can be located at network discontinuities only (merges, diverges, and lane-drops). In this way, the spatial discontinuities can address easily. The computation cost depends on the number of cell boundaries (x-dimension) in the network and the number of vehicles (n-dimension) to propagate during the simulation. Therefore, this would improve computational efficiency for large scale applications.

More importantly, this mesoscopic scheme is particularly convenient for data assimilation.

In reality, the flow characteristics are mostly observed at fixed points (e.g., spatial fixed loop data)

or along vehicle trajectories (e.g., vehicle-number fixed probe data). As discussed in literature that the Eulerian formulation is suitable for incorporating loop data and the Lagrangian-time formulation is suitable for probe data assimilation, the Lagrangian-space formulation is considered to be well-compatible for assimilating both types of observations. Because these observations are located on cell boundaries of the mesoscopic grid, which makes any traffic state estimation method convenient with this approach/formulation. This formulation can be easily coupled with any data assimilation techniques to perform state estimation. Due to the nature of the mesoscopic system model, the TSE might be not restricted to discretized mesoscopic x -n grids. If we know any two boundaries in the network and an observation at a certain location or of a certain vehicle, we can generalize TSE for this specific assimilation problem.

METHODOLOGY

Methods for estimating traffic states based on loop and probe data are presented in sections 4.1 and 4.2, respectively. Next a method that combines both data sources is presented in section 4.3.

First, three definitions with respect to different traffic states are given in the following: As mentioned by the authors in (3), the update of the model is a parsimonious adjustment of the demand and/or the supply terms at cell boundaries. It has to be implemented so that the CFL stability condition is respected. The reader is referred to the paper for more details and validation.

• an observation (o-)

TSE based on probe data

The data assimilation framework presented above is limited to Eulerian (loop) data while nowadays increasing amount of traffic data are collected by Lagrangian (probe) sensors. Thus a TSE estimator based on Lagrangian observations becomes essential for real applications. Probe sensors collect positions of equipped vehicles at a given time frequency. They are usually processed for providing aggregated indicators, for instance the mean speed per link. However, most of the wealth of probe data is lost during the aggregating process.

In this paper, the TSE estimator enables to assimilate positions and times without any aggregation process, which allow for using most of the details of probe data. The method is divided into 4 steps.

• and N b (t, x down ) can then be considered for estimating n-index of the probe based on variational principles, as illustrated in Figure 2.

n a p,i = min N b up,u , N b down,w + k x .(x down -x o p,i ) (9)
where

     N b up,u = N b t o p,i - x up u , x up N b down,w = N b t o p,i - x down -x o p,i w , x down
Equation 9 provides the n-index estimated locally (for a single time-position). At this stage, local n-index estimation could be flawed by four sources of errors: errors on the model parameters, errors on the boundary conditions, non-FIFO traffic conditions or occurrence of a traffic incident. Local errors on the estimated n-index may induce global inconsistencies on the resulting arrival/supply times. To tackle this problem, a global optimization is developed and it consists of two steps. The first step aims at building the variational proximity matrix , which returns the variational cost (in veh.) between each of the time-space observations from probes (with respect to the variational principles (1)). Based on that, the second step calculates the optimal n-index, denoted n * . The optimal solution minimizes the entropy of the system while keeping a constant n-index along probe trajectories. The optimization procedure searches in the range of all possible n-indices, and this search range is defined by the minimum and maximum values from the variational principles and the range of local n-index estimation. The entropy is defined as:

E(n * p ) = i n * p,i n local p,i .ln( n * p,i n local p,i ) (10)
The final solution consists of the triplets {n * p , x o p , t o p }, where n * p is the optimal n-index, and

t a a,n a p = t o p,i + x down -x o p,i u (11)
Upstream: update of supply times The upstream cell boundary is influenced by probe vehicles located in a time window with a length P and that moves with a maximum jam speed w, see Figure 3(c). For each probe vehicle, triplets {n a p,i , x o p,i , t o p,i } are considered as internal boundary conditions to revise supply times at the cell boundary upstream. Within the influencing area, the updated supply times respect as illustrated in Figure 3(d):

t a s,n a p +(x o p,i -xup).kx = t o p,i + x p p,i -x up w , ∀i ∈ S p (12)
CFL condition The data assimilation process is sequential with time steps based on data time frequency ∆T . The CFL stability condition has to be respected during the sequential update of the traffic model. It requires that each cell boundary has to be updated over a time period ∆T U , which is bounded as a wave cannot travel through a whole cell during this time period. Consequently, if ∆T ≥ ∆T U then the updating process must proceed step by step (as described in the previous section) with a maximum time step ∆T U .

Assimilating both loop and probe data

Loop and probe data provide information of different nature, it is therefore impossible to fuse the two data sources to perform a one-shot assimilation process. Reviewing their respective actions, the two TSE estimators act in a complementary manner. On one hand, TSE based on loop data allows for an adjustment of the flow by adding -deleting -advancing -delaying vehicles at loop On the other hand, the TSE based on probe data adjusts arrival and supply times at cell boundaries considering probe trajectories as internal cell-boundary conditions. From a physical point of view, it acts as a 'travel time regulator' along cells travelled by probe vehicles.

To make the best potential use of both data, we propose first to estimate trafic states from loop data at loop sensors locations and then to estimate traffic states from probe data everywhere else. The main reason for this sequence is the following : TSE based on loop data improves the flow estimation at cell boundaries and therefore enhance the TSE based on probe data along cells.

It results in a 7 steps methodology:

•

MODEL VALIDATION

This section aims to analyze and validate the performance of the TSE methodology with loop sensors and probe sensors (separately and jointly).

Experimental validation setup

The ground truth is emulated based on a microscopic LWR model (Newell's car-following model, equivalent to the LWR model at a macroscopic scale). The model has been run on a homogeneous road stretch (L = 2000m, single lane) with a demand-supply scenario so that a congestion propagates through the network, see vehicle trajectories in Figure 3. A loop sensor located in the middle of the network (x = 1000m) collects flows and speeds with an aggregation period of 1-minute.

Moreover, 5% of the vehicles are considered as probe sensors for which time-position information is reported at every 30s.

The traffic flow model is a LS-LWR model. The network is composed of 2 cells of 1000m in length, upstream and downstream of the loop sensor location. The demand-supply scenario has also been predefined with an approximative demand and a high supply so that traffic conditions are always free-flowing on the network. In summary, when a congestion occurs, loop sensors can estimate travel times providing that congestion states have passed over the loops. The result shows that travel times might be underestimated over the network level. And this underestimation will become significant when traffic congestion is triggered far downstream the loop sensor. We conclude that for operational purposes loop sensors have to be located as close as the triggering location of a jam/bottleneck to provide accurate estimation. In addition, the complementary information from downstream loop sensors can improve the performance of data assimilation.

TSE based on loop observation model

TSE based on probe observation model

Figure 5 provides the estimated traffic states considering probe data only. The performance of TSE based on solo probe data provides similar performance over the two cells. It is noteworthy that TSE is very responsive as the congestion phenomenon occurs, mainly due to the probe data with a homogeneous coverage of the network both in time and space. It should also be noted that travel times are underestimated in this validation scenario due to the experimental setup. The traffic model considers a low demand versus high supply scenario. Information from probes allows for an adjustment of the supply times at the intercell boundary, but it does not rectify the underestimated flow demand (from downstream) and thus underestimate travel times. Note that the result We conclude that for operational purposes, the knowledge of the demand at any point of the network is decisive and critical when probe data are used for estimating traffic states. The estimation can be enhanced with an accurate prior estimation of the demand; or combining probe data with loop data, as proposed in section 4.3.

TSE based on loop and probe observation model

Downstream cell

The performance is slightly enhanced compared to the results obtained with probe observations only. It confirms that both observations are very complementary when assimilated in the framework proposed in the paper.

APPLICATION TO A FREEWAY CORRIDOR

The previous section demonstrates the exactness of the estimator when applying to a network with FIFO conditions and homogeneous driving behavior. These assumptions are restrictive and not reflective of reality. This section aims at evaluating the performance of the estimator considering a multi-lane corridor with on-and off-ramps, with a relaxed FIFO assumption and distributed driving behavior.

Preparation of the observational model

Ground truth data have been emulated based on a microscopic traffic simulator (FOSIM (2)). This simulator is developed at the Delft University of Technology, specially designed for the detailed analysis in freeway networks. All the parameters in terms of driving behaviors have been calibrated and validated based on data from Dutch freeways. A three-lane freeway with one on-ramp and one off-ramp is designed, as illustrated in Figure 7(a) (the first 500 m as the warming-up section in Simulation, the last 1000m as the cooling-down section).

A demand-supply scenario has been built in such a way that a congestion is onset at the onramp. The model has been run twice: scenario 1 provides traffic conditions with only passenger cars whereas scenario 2 considers a mixed traffic condition (with 90% cars, 10% trucks). The resulting time-space diagrams and travel times are illustrated in Figure 7.

Based on FOSIM simulation results, Eulerian and Lagrangian observation models have been built. First, loop sensors have been located on the main road : loop 1 -100m after the entrance of the network, and loop 2 -100m upstream of on-ramp. Second, 10% equipped probe vehicles return their exact positions every 20s.

Preparation of the traffic model

The traffic model is the mesoscopic LWR model applied on a network with 7 cells : five cells for the main road (numbered from 1 to 5), one cell 6 for the off-ramp and one cell 7 for the on-ramp. Cell As excepted, TSE combining loop and probe observations outperforms the estimations considering loop and probe separately. The increase of travel times is detected immediately after the onset of congestion and the estimated travel times dynamically correspond to the ground truth. sensors are spatially fixed, they only catch supply information as information propagates upstream.

Hence, travel time estimation might be inconsistent (see section 5.2, the case of the downstream cell). It is therefore critical to locate loop sensors at the spots of a jam and a bottleneck to provide accurate estimation.

In contrast, Lagrangian observations spread over the network in space and time. Given a reasonable resolution (10% in the validation cases), they can update both arrival times and supply times without any latency. However, two limitations exist in this method. First, the n-index estimation relies on the assumption of the FIFO condition, which is unrealistic. This might lead to poor estimation of travel times when probe information deviates from the average traffic conditions (see section 6.3). Second, the demand/flow at cell boundaries cannot be adjusted. A prior estimation of the demand will improve the performance of TSE with probe data (see section 5.2).

The combination of the two data sources compensates the limitations of each other. The experiment results demonstrate that TSE with data combination outperforms the estimation with a single source. Lagrangian observations provide a homogeneous coverage of the network in time and space for a discrete set of (probe) vehicles. Probe vehicles allow for an revision of demand-supply times at neighboring cell boundary of the network. The update acts as a cell 'travel time regulator' that yields good results under the condition that the demand on the network is known. Note that the nindex values of probe vehicles are critical and essential in the proposed approach. The calculation of this variable is application-specific (under FIFO or non-FIFO condition) regarding estimation performance, and it is subjective for further investigation.

CONCLUSION

In the TSE framework with data combination as proposed in section 4.

FIGURE 1

 1 FIGURE 1 Numerical solutions in Lagrangian-space coordinates

  Step 1: the o-state and the b-state are collected and transformed • Step 2: Global Analysis, which consists of estimating the n-index of probe vehicles • Step 3: the model is updated accordingly, which consists of adjusting arrival and supply times at cell boundaries of the model • Step 4: the model is run over the next sequence to provide a new background The two following sections elaborate steps 2 and 3, which are the keys to successfully update traffic states.

FIGURE 2

 2 FIGURE 2 n-index estimation

t o p and x o p are theFocus on step 3 :

 3 observed time and position of the probe p. Update of arrival and supply times at cell boundaries Once a-states are known, probe trajectories are considered as internal cell boundary conditions that are transformed into demand or supply conditions at neighboring cell boundaries. Here, we present the update of the arrival and the supply times at a cell boundary over a period P , considering that a set of probe vehicles has been analyzed.Downstream: update of arrival timesThe downstream cell boundary is influenced by probe vehicles located in a time window with a length P and that moves with a free-flowing wave speed u, see Figure3(a). Within the influencing area, each probe vehicle provides information on its upcoming arrival times. When probe vehicles travel through a cell, successive time-positions provide feasible arrival times at the downstream cell boundary. For each probe vehicle, only the latest triplet {n a p , x o p,i , t o p,i } is considered for updating the arrival time at the cell downstream, as illustrated in Figure3(b).

  Update of the supply time sensors locations. From a physical point of view, it acts as a 'flow regulator' at cell boundaries.

  Step 1: collection and transformation of the loop data and the model background states • Step 2: Global Analysis, which consists of estimating headway-regime pairs (a-states) at each loop location • Step 3: the model is updated accordingly (see section 4.1). At this stage, the updated model provides the best possible estimated traffic states at cell boundaries. This version of the model is considered as a new model background to be combined with probe data • Step 4: collection and transformation of the probe data and the (updated) model background • Step 5: Global Analysis, which consists of estimation the n-index of probes along cells • Step 6: update of the model accordingly, by revising arrival and supply times at every cell boundaries, except those already updated during the step 3. • Step 7: run the model over the next sequence Here again, this sequence has to be implemented while respecting CFL stability condition mentioned in section 4.2.2.

Figure 4

 4 Figure 4 provides the estimated traffic states considering data from the loop sensor. In this figure, traffic states have been rearranged to provide travel times over the two cells. The red line provides the reference (ground truth) travel times and the blue line returns the reconstructed travel times.

FIGURE 3

 3 FIGURE 3 Observational model

FIGURE 4 FIGURE 5

 45 FIGURE 4 Travel time estimation from loop data assimilation

Figure 6

 6 Figure 6 provides the estimated traffic states considering loop and probe data. The results show the travel times estimated here over performs and cumulates the benefits mentioned for loop and probe observation model considered separately.

FIGURE 6

 6 FIGURE 6 Travel time estimation from loop and probe data assimilation

  time-space diagram (c) time-space diagram (d) travel times (e) travel times

FIGURE 7 FOSIM

 7 FIGURE 7 FOSIM observation models: cars only (a and c) and mixed traffic (b and d)

FIGURE 8

 8 FIGURE 8 Comparison of estimated travel times

  Main findingsA TSE estimator based Eulerian observations combined with a mesoscopic LWR model has been proposed and validated in (3). This paper complements the methodology with Lagrangian observations. Now both Eulerian and Lagrangian observations can be used for TSE in a unique framework.Eulerian observations provide comprehensive observations in time and vehicle for a discrete set of locations in the network. At those locations, model states are successfully revised, which provide good performance when observations are located near the head of congestion. The update acts as a 'flow regulator' at cell boundaries by adding, deleting, advancing or delaying vehicles.

3 ,

 3 Eulerian and La-grangian observations become highly complementary. Eulerian observations successfully update traffic states (especially the flow) at loop locations of the network while Lagrangian observations successfully update cell travel times along the network. The methodology has been verified on the synthetic data derived from the same underlying traffic flow model. The proposed TSE framework has been applied to a freeway corridor with a relaxed FIFO condition and distributed driving behavior. The validity and performance have been tested using the ground truth from a microscopic simulator. Further research Future research includes, (a) to test the framework in a more general case, a realistic large-scale network with multiclass traffic, and a non-FIFO condition; (b) to test different sources of data observations, e.g., bluetooth data, information from connected vehicles; (c) to apply specific data assimilation techniques to account for model and observation reliability; (d) to study the optimal layout of Eulerian sensors and the minimal penetration of Lagrangian sensors for accurate estimation; (e) online estimation of traffic demand and important parameters in the traffic flow model; (f) traffic state prediction based on the proposed estimation framework.

  A data assimilation method using sole loop data first proposed in (3). It requires the numerical scheme to be set as follows: ∆n to 1 and cells boundaries at each loop location. It considers flow and speed time series collected by loop sensors at locations {X loop } with a given frequency ∆T . Then it is implemented as a sequential procedure, for which each sequence is divided into 4

	state is a traffic state measured by a sensor
	• a background (b-) state is a state forecasted by a traffic flow model
	• an analysis (a-) state is the result of an analysis procedure (or algorithm) that provides
	the most likely state regarding o-and b-states
	TSE based on loop data
	successive steps:

• Step 1: the o-state and b-state are collected and transformed • Step 2: a Global Analysis is performed to (a-) state • Step 3: the state of the model is updated accordingly, by adjusting arrival and supply times at cells boundaries • Step 4: the model is run to provide a background state for the next sequence

Table 1

 1

	provides three Measurements of Effectiveness (MoEs) that have been calculated
	over the period t = [20 -45]min: mean average error (MAE) , mean average percentage error
	(MAPE) and mean percentage error (MPE). All the MoEs globally confirm the previous comments.

DISCUSSION

Based on the previous results, we conclude that Eulerian observations can update arrival and supply times, by adjusting the demand (flow) via adding or deleting vehicles locally. However, since loop

TABLE 1

 1 Performance of the different observational models over the indicated simulation period, Scenario 1 (cars only) And Scenario 2 (cars and trucks)

	SCENARIO 1 model model + loop model + probe model + (loop and probe)
	RMSE (s)	96	69	28	21
	MAPE (%)	32	24	9	6
	MPE (%)	-32	-24	1.43	1
	SCENARIO 2 model model + loop model + probe model + (loop and probe)
	RMSE (s)	229	92	56	42
	MAPE (%)	49	22	12	10
	MPE (%)	-49	-22	4	3

ACKNOWLEDGMENT

The authors sincerely thank L. Leclercq (University of Lyon, ENTPE, IFSTTAR, LICIT, UMR-T9401) for the fruitful discussions and assistance on the methodology presented in the paper.