Practical Model Predictive Control for a Class of Nonlinear Systems Using Linear Parameter-Varying Representations

(1)

Practical Model Predictive Control for a Class of Nonlinear Systems Using Linear Parameter-Varying Representations

HOSSAM S. ABBAS ^1,2, (Senior Member, IEEE), PABLO S. G. CISNEROS ³, GEORG MÄNNEL¹, (Graduate Student Member, IEEE), PHILIPP ROSTALSKI ¹, AND HERBERT WERNER³

1Institute for Electrical Engineering in Medicine, University of Lübeck, 23558 Lübeck, Germany 2Electrical Engineering Department, Faculty of Engineering, Assiut University, Assiut 71515, Egypt 3Institute of Control Systems, Hamburg University of Technology, 21073 Hamburg, Germany

Corresponding author: Hossam S. Abbas (h.abbas@uni-luebeck.de)

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Project 419290163.

ABSTRACT In this paper, a practical model predictive control (MPC) for tracking desired reference trajectories is demonstrated for controlling a class of nonlinear systems subject to constraints, which comprises diverse mechanical applications. Owing to thelinear parameter-varying(LPV) formulation of the associated nonlinear dynamics, the online MPC optimization problem is solvable as a singlequadratic programming(QP) problem of complexity similar to that of LTI systems. For offset-free tracking, based on the notion ofadmissible reference, the controller ensures convergence to any admissible reference while its deviation from the desired reference is penalized in the stage cost of the optimization problem. This mechanism provides a safety feature under the physical limitations of the system. To guarantee stability and recursive feasibility, a terminal cost as a tracking error penalty term and a terminal constraint associated with both the terminal state and the admissible reference are included. We use tube-based concept to deal with the uncertainty of the scheduling parameter over the prediction horizon. Therefore, the online optimization problem is solved for only the nominal system corresponding to the current value of the scheduling parameter and subject to tightened constraint sets. The proposed approach has been implemented successfully in real-time onto a robotic manipulator, the experimental results illustrates its efficiency and practicality.

INDEX TERMS Constrained systems, linear parameter-varying systems, model predictive control, robust stability, robotic manipulators.

I. INTRODUCTION

The ultimate goal of a control system is to achieve stability and a desired level of performance for plants which often have nonlinear(NL) dynamics, constrained levels of operation and are subjected to disturbances and measurement noise.Model predictive control(MPC) [1] is a paradigm that can system- atically handle such complications. It is a control approach which can optimize system performance online based on predicting its future behavior over a so-called prediction horizon(N). However, unless the MPC optimization problem is formulated appropriately, computational complexity or inherent conservatism could affect performance or even result in infeasibility or instability.

The associate editor coordinating the review of this manuscript and approving it for publication was Zheng Chen .

Linear parameter-varying(LPV) approach [2] is a promis- ing framework for controlling nonlinear and time-varying (TV) systems using linear control techniques. Numerous successful industrial applications have proven its efficacy, see, e.g., [3], [4]. Adiscrete-time LPV system is represented in state-space form as

x(k+1)=A(p(k))x(k)+B(p(k))u(k), (1a) y(k)=C(p(k))x(k)+D(p(k))u(k), (1b) whereu(k) ∈ Rⁿ^u,x(k) ∈ Rⁿ^x, y(k) ∈ Rⁿ^y andp(k) ∈ Rⁿ^p, are vectors of the input, state, output and scheduling variable (parameter) of the system at a time indexk∈Nand A,B,C,Dare parameter-dependent matrices with appropriate dimensions. The representation in (1) provides a modeling framework that can efficiently describe NL/TV systems in a

(2)

linear setting, in which the relation between the input and output signals is linear, but dependent on p. It is assumed that the scheduling variable is measurable and taking values from a so-called scheduling range P. In LPV models of NL applications, p is often associated with the input, state or output (endogenous signals) of the system, thus, LPV representations of NL systems are often referred to asquasi- LPV (qLPV) models. The LPV system representation (1) is usually employed for controller design based on linear optimal and robust control methods [5]. In all LPV control strategies, closed-loop stability and performance guarantees during implementation are established under the assumption thatpwill stay all the time insideP. However, in case of qLPV representations, wherepis not a free variable, i.e., endogenous, such assumption cannot be ensured, unless the control design methodology can restrictpfrom deviation outsideP via input/state constraints. In fact, most of the LPV control strategies based on H₂/H∞ optimal control, e.g., [6], [7], cannot handle signal constraints, which renders achieved stability/performance guarantees based on these strategies void ifpviolates its region. On the contrary, MPC can optimally tackle such a difficulty by virtue of constraints handling.

Although generalized formulations of MPC for NL systems exist, in the context of nonlinear MPC (NMPC), MPC schemes based on LPV systems (LPVMPC) for controlling NL plants have become popular in the control community since the late 1990s, see [8]–[11] and the recent survey in [12]. LPVMPC provides an intermediate step between conventional linear MPC (LMPC) and NMPC in terms of the achievable control performance with a moderate computational complexity. Usually the achievable control performance of an MPC design approach reflects its degree of conservatism. The strength of the LPVMPC is its ability to be solved using LMPC tools while it can achieve asymp- totic stability andrecursive feasibility¹ guarantees for NL systems, yet with a computational cost that is lower than NMPC [11].

However, the major difficulty of an LPVMPC setting is that the scheduling variable is accessible only at the current time instantk, but its future values, which are required over the prediction horizon of the MPC for state prediction, are unknown. Therefore, robust MPC is used to handle such uncertainty of p; however, it is usually conservative ifp is assumed to vary arbitrarily fast over its full range, see, e.g., [13] and [14], which might affect achievable performance.

Recently, tube-based MPC [15] has been investigated for LPV systems in [16], [17] and [18] to handle the uncertainty of the future values of p by exploiting known bounds on its rate of variation to obtain an admissible range of pover N instead of considering its full range, which can lead to low conservative techniques. However, the computational complexity of the associated optimization problem can be very costly for practical applications. Note that one of the

1In MPC, recursive feasibility is defined as follows: If the MPC optimization problem is initially feasible, then it will remain feasible.

main challenges of any MPC algorithm is to compromise the degree of conservatism with the computational complexity.

In this paper, to tackle these difficulties, a tube-based setting is carried out, where the bounds on the rate of variation ofpcan be exploited to construct scheduling tubes containing the N possible future values of p, which can considerably reduce the conservatism of considering its full range. Such tubes are employed to construct state tubes as rigid tubes [15]

to which the future trajectories of the state are confined and used for constraint tightening. In contrast to many LPVMPC approaches in the literature, the computation of these tubes in our strategy is performed offline, which significantly reduces the online computations.

To put any LPVMPC algorithm for practical use, it should be suited for reference tracking. Such a control problem has been rarely investigated in the context of LPVMPC, where most of the developed methods have focused on regulating the state of the controlled system to the origin or to a set point. Robust tracking MPC approaches for LPV systems have been developed recently in [19], which have achieved offset-free tracking for piecewise constant references. In [20]–[22] LPVMPC algorithms have been proposed for tracking time-varying reference trajectories, e.g., command trajectories in robotics applications. However, all these approaches cannot guarantee recursive feasibility of the associated MPC optimization problem.

For practical implementation, we extend our proposed MPC strategy to include reference tracking. Inspired by the approach of [23] and [24], for a given desired reference trajectory, the corresponding admissible steady sate and input are parameterized by a parameter vector referred to as the admissible output, which is among the decision variables of the optimization problem, and its deviation from the desired reference is penalized in the MPC cost function. This can lead to offset-free tracking if the desired reference is admissible, then, the system is steered toward the closest admissible reference. Such mechanism allows a safety constraint, which protects the system from tracking references beyond its physical limits without external interruption of its operation.

Moreover, a larger domain of attraction can be achieved in comparison with standard MPC for tracking. Moreover, at the expense of increasing the number of decision variables associated with the admissible output, the proposed approach is not restricted for tracking piecewise constant references as that of [23]. For guaranteeing recursive feasibility, we utilize the concept of invariant set for tracking which is associated with the maximum set of admissible references. That invariant set is employed as a terminal set in the LMPVMPC problem for reference tracking, thus, asymptotic stability of the closed-loop system can be ensured.

To validate the practicality of the proposed technique, it is implemented experimentally on a twodegree-of-freedom (DOF) robotic manipulator for reference tracking. The unique property of the proposed techniques is that it solves the online LPVMPC optimization problem for the nominal system corresponding to the current value ofp. At the expense of losing

(3)

some performance, such formulation can result in computational burden comparable to that of conventional LMPC.

To summarize, the contributions of this paper are as follows:

1) A general LPV modeling which can accommodate a wide class of mechanical systems that can be represented by rigid body nonlinear dynamics. This allows a straight forward formulation for the MPC reference tracking problem.

2) A computationally convenient LPVMPC algorithm based on quadratic programming (QP), which can ensure offset-free tracking and is appropriate for controlling several nonlinear applications.

3) Guarantees for recursive feasibility of the LPVMPC optimization problem and asymptotic stability of the closed-loop system.

4) A successful experimental implementation to control a 2DOF robotic manipulator for reference tracking, which shows comparable results to a recently developed more computationally demanding MPC scheme.

The paper is organized as follows: To introduce the proposed approach, some preliminaries from [17] are given in SectionII. The modeling aspects related to the applications considered in this work are illustrated in SectionIII. The proposed LPVMPC approach for reference tracking is developed in SectionIVtogether with other related computation issues.

The application to the 2DOF robotic manipulator and the experimental results are demonstrated in SectionV. Finally, in SectionVIthe conclusion is given.

A. NOTATION AND DEFINITION

LetNdenote the set of non-negative integers including zero.

We denote the predicted values of a variablex(k) at timek+i based on the available information at timekasxi|ksuch that x0|k = x(k). Co{·}denotes the convex hull of a set. For any vector x ∈ Rⁿ,kxk denotes the 2-norm, and the weighted norm is defined bykxk²

P=x^>Px, whereP=P^>,P∈R^n×n. A polytope is a compact polyhedron, which is the inter- section of a finite number of half spaces. A (hyper)box is a convex polytope where all the defining hyperplanes are axis parallel. A PC-set is a set that is convex, compact and has a nonempty interior containing the origin. Given two sets A⊂RⁿandB⊂Rⁿ, the Minkowski set addition is defined byA⊕B := {a+b | a ∈ A,b ∈ B}and the Pontryagin set difference is defined byA B:= {a|a⊕B⊆A}. Let a∈ Rⁿâ,b∈Rⁿ^b and a setϒ ⊂Rⁿâ⁺ⁿ^b, then the projection of ϒ ontoa is defined asProj_a(ϒ) = {a ∈ Rⁿâ | ∃ b ∈ Rⁿ^b,(a,b)∈ϒ}.

A functionf :R⁺→R⁺is of classK∞if it is continuous, strictly increasing,f(0)=0 and lim_ξ→∞f(ξ)= ∞.

Definition 1 (Quadratic Stabilizability [25]): The system (1) is quadratically stabilizable if there exists a positive definite function V : x → x^>Px, where P = P^> 0, P ∈ Rⁿ^x^×n^x and a control law u = Kx,K ∈ Rⁿ^u^×n^x such that

V A^c(p)x

−V(x)≤ −x^>(Q+K^>RK)x (2)

∀x∈Rⁿ^x,p∈P, whereQ=Q^>0,R=R^>0 and A^c(p)=A(p)+B(p)K. (3) Then the origin is globally exponentially stable forx(k+1)= A^c(p(k))x(k),∀p∈P.

Definition 2 (Robust Positive Invariant Set [26]): For the system (1), with state and input constraint sets X and U, respectively, and the control law u(k) = Kx(k) ∈ U, the setXf ⊂ Xis robustly positively invariant (RPI) if for all x(k)∈Xfandp(k)∈P,x(k+1)∈Xf.

Definition 3 [Robust Invariant Set for Tracking (RIST)]:

It is the set of all initial states and steady states and inputs of the system (1) that can be stabilized by the control law

u(k)=K(x(k)− ¯x)+ ¯u, (4) wherex¯andu¯denote the steady state and input, respectively, andK yields the closed-loop system matrixA^c(p) as shown in (3) Schur for all values ofp ∈ P. Moreover, the control law (4) fulfills the input constraints and renders the system state constraints satisfied throughout its evolution.

Definition3is extended from [24].

II. PRELIMINARIES

In this section, we review some material from [17]. Then, we introduce a low complexity MPC scheme for regulating the state of an LPV controlled system into the origin, which will be extend in the sequel of the paper into an MPC for tracking a given reference trajectory. The latter will be implemented for controlling the robotic manipulator application in this work, which can be applicable also for a broad class of mechanical systems.

Consider discrete-time LPV systems represented by (1) and let the following assumptions be satisfied.

Assumption 4: (i) The system (1) is quadratically stabilizable.

(ii) The values ofx(k) andp(k) are available at every time k∈N.

(iii) The setsXandUare polytopic PC-sets.

(iv) The system matrices depend affinely onp, i.e., A(p)=A⁰+

np

X

j=1

p^jA^j, B(p)=B⁰+

np

X

j=1

p^jB^j, (5) wherep^jdenotes thej^th-element of the vectorpandA^j, B^j, are constant matrices with appropriate dimensions.

(v) The parameter setPis a compact hyper-box defined as P:= {p∈Rⁿ^p | p^j^,^min≤p^j≤p^j^,^max, j=1,· · ·,np}. (vi) The rate of variation ofpis denoted as

dp(k)=p(k)−p(k−1), which is bounded such that dp∈dP, where

dP:= {dp∈Rⁿ^p | |dp^j| ≤dp^j^,^max, j=1,· · · ,np} is a compact hyper-box.

(4)

The above assumptions are imposed for the mathematical derivations; however, they are standard in LPV literature and encompass many practical situations [6].

In order to quantify the uncertainty ofpover the prediction horizon, letpi|k denote the uncertain scheduling parameter at any step iover N such that pi|k ∈ P_i|k, see Fig.1, where P_i|k ⊆Pis a compact hyper-box defined as

P_i|k :=

(

p_i|k∈Rⁿ^p | |p^j_i|k− ˘p^j_i|k| ≤

i

X

l=1

dp^j_l|k^,^max,

j=1,· · ·,n_p

. where dp^j_l|k^,^max≤dp^j,^maxandp˘i|k is a known (nominal) value of the scheduling parameter, which can be computed at the center ofP_i|k, see Fig.1for illustration. Thereforep_i|kcan be parameterized as

p_i|k= ˘p_i|k+

i

X

l=1

dp_l|k, (6)

where dp_l|k∈dP.

Define the uncertain system over the prediction horizon as xi+1|k =A(pi|k)xi|k+B(pi|k)ui|k, pi|k ∈P_i|k, (7) for all i = 1,· · ·,N − 1, where x_i|k represents the corresponding uncertain state at stepi. We can rewrite the system (7) as

xi+1|k =Ai|kxi|k+Bi|kui|k+wi|k, (8) whereA_i|k =A(p˘_i|k),B_i|k =B(p˘_i|k) are known matrices and

w_i|k :=(A(p_i|k)−A_i|k)x_i|k+(B(p_i|k)−B_i|k)u_i|k, wherew_i|k represents an additive disturbance describing the uncertainty ofpi|k such thatwi|k ∈Wi|k,

W_i|k :=Co

A(p_i|k)−A_i|k

x_i|k+ B(p_i|k)−B_i|k u_i|k,

| pi|k ∈P_i|k, xi|k ∈X, ui|k ∈U , (9) for all i = 1,· · · ,N −1. Any set W_i|k is polytopic and contains the origin, see [17]. Based on this setting the value ofx1|k in (8) at any timek ≥ 0 is known asx0|k = x(k) is assumed to be known, thus,w0|k =0.

Introduce the LPV nominal system for (8) as

zi+1|k =Ai|kzi|k+Bi|kvi|k. (10) wherez_i|k represents the nominal state at stepi,v_i|k is the nominal input, which is related toui|k in (8) by

u_i|k =v_i|k+K(x_i|k−z_i|k), (11) whereK ∈ Rⁿ^u^×n^x is referred to as the disturbance con- troller[15], which is used to penalize the error betweenx_i|k andz_i|k. For simplicity, we consider a robust controller K according to Definition1. Letz0|k =x0|k, thus,z1|k =x1|k

asw0|k =0, compare (8) and (10) wheni=0. Furthermore, introduceεi|k =xi|k−zi|k to denote the difference between the uncertain and the nominal states with the dynamics

εi+1|k =A^c_i_|_kεi|k+wi|k, (12) whereA^c_i|k =A_i|k +B_i|kK is a Schur matrix for anyi,kand wi|k ∈W_i|k. Moreover, letεi|k ∈S_i|k, where

S_i|k:=W_i−1|k⊕A^c_i−1|kS_i−1|k, i=2,· · · ,N, (13) which is a PC-set [17]; note thatε0|k =ε1|k =0 and hence S0|k=S1|k = {0}.

The above formulation implies that the state trajectories overNare confined in a state tube with varying cross section and shape according toS_i|k, which sometimes is referred to as heterogeneoustube [18]. Similarly, every control trajectory over N lies in a control tube with cross sections KS_i|k. Constructing heterogeneous tubes yields less conservative tubes in comparison with homothetic or rigid tubes [18];

however, that is at the expense of considerable increase the MPC computational burden.

Based on the above formulation, a tube-based MPC for LPV systems has been introduced in [17]. Given the initial conditionsz_0|k =x_0|k = x(k) andp_0|k =p(k), the optimal values of the nominal statez_i|k, for alli = 1,· · ·,N −1, and the nominal control inputv_i|k, for alli=0,· · ·,N−1, can be computed at anyk ∈ N, such that the state and input constraints of the LPV system (8) are satisfied by solving the following optimization problem

v0|k,···min,vN−1|k

N−1

X

i=0

kz_i|kk²

Q+ kv_i|kk²

z_N|k ⊂Xf S_N|k (14d)

and the nominal system dynamics (10) with z_0|k = x(k), where Q = Q^> 0 and R = R^> 0 are used as tuning parameters to meet some desired performance in the stage cost of (14a) andV_f(·) is the terminal cost, which is chosen offline as well as the terminal setXf ⊂ Xin (14) to guarantee asymptotic stability of the closed-loop system and recursive feasibility of the above optimization problem.

The computations ofVf(·) together with a terminal controller K andXf as an RPI set underK are carried out according

(5)

to Definitions 1 and 2, respectively. Therefore, according to the basic idea of tube-based MPC [15], the satisfaction of the constraints (14b-d) ensures thatxi|k ∈ Xfor alli = 1,· · · ,N −1, u_i|k ∈ Ufor all i = 0,1,· · · ,N −1 and x_N|k ∈Xf, see [17] for more details.

The online implementation of (14) at each timekincludes also computingW_i|k for alli=1,· · · ,N−1 using (9) and S_i|k for all i = 2,· · ·,N using (13) as well as performing the constraints tightening in (14b-d). Then, the optimization problem (14) can solved as a QP problem. However, the online computation burden of computing W_i|k andS_i|k might be relatively higher than solving the MPC optimization problem due to the involved sets addition/subtraction.

For practical implementation, we modify in the following the above MPC approach of [17] —at the expense of a prob- able increase of conservatism— to avoid the online computations of the setsW_i|k,S_i|kand the associated inline tightening of the constraint sets. Consider the following assumption:

Assumption 5: There exists a PC-setS⊂Xfsuch that S_N|k ⊆S (15) holds for allk∈N.

Note that the conditionS⊂Xfin Assumption5is necessary forZfto have an interior. We will discuss in SectionIV-D the computation ofS(offline), which can satisfy (15). Based on Assumption5, the tightened constraint sets in (14b-d) can be replaced by the following sets

z_1|k ⊂X, (16a)

zi|k ⊂Z, i=2,· · ·,N−1, (16b)

v0|k,v1|k ⊂U, (16c)

v_i|k ⊂V, i=2,· · ·,N−1, (16d)

zN|k ⊂Zf, (16e)

where

Z⊂ X S, (17a)

V⊂ U KS, (17b)

Zf⊂ Xf S, Zf⊂Z. (17c)

Now, given the setS, all tightened constraint setsZ,VandZf

can be computed offline, thus, we avoid the online computation of the setsW_i|kandS_i|k.

Remark 6: Such formulation implies that the state trajectories over the prediction horizon are confined in a state tube with a center z_i|k and a fixed cross section S, i.e., x_i|k ∈ z_i|k⊕S,∀i=0,1,· · · ,N, which sometimes is referred to as rigid tube [27]. Similarly, every control trajectory overNlies in a control tube asu_i|k ∈v_i|k ⊕KS. Note that the feedback policy in (11) can affect the size ofS.

Therefore, the online computation involves just updating the system matrices based on the current value of the scheduling variablep(k) and solving the optimization problem

v_0|k,···min,v_N_−1|k N−1

X

i=0

kz_i|kk²_Q+ kv_i|kk²_R+Vf(z_N|k) (18a)

subject to (16a-e). (18b)

To construct the stage cost in (18a) fori=0,1,· · ·,N−1, one can consider the LPV nominal model (10) or just a simplified LTI model given by

z_i+1|k =A_0|kz_i|k+B_0|kv_i|k, (19) for alli = 0,1,· · · ,N −1, withz_0|k = x(k), in this case, the uncertainty ofp_i|k is parameterized as (6) wherep˘_i|k is replaced byp0|k. Considering (10) might give better information about the evolution ofpover the prediction horizon;

however, it requires the computation of the nominal scheduling variablep˘i|k overN, which takes into account the upper bound on the rate of change ofp. The proposed simplified version of the MPC compared with that of [17] can result in a lower computational complexity, which is comparable to that of conventional MPC for LTI systems; however, depending on the size of the setS, this might be more conservative than that in [17].

Moreover, at the end of the prediction horizon, i.e., at step N, the state should satisfy the terminal constraintx_N|k ∈Xf. This condition is necessary for guaranteeing stability and recursive feasibility. Satisfaction ofx_N|k ∈ Xfis ensured if the nominal system at stepN satisfies the tightened terminal constraint, i.e.,z_N|k ∈ Zf, this holds true due to the condition (17c), wherex_N|k ∈ z_N|k ⊕S. Sincex_i|k ∈ Xf for all i≥N ifxN|k ∈Xf—asXfis RPI for (1) under the controller K— we can conclude thatzi|k ∈Zffor alli≥NifzN|k ∈Zf. Moreover, sinceKXf ⊂Uholds for allp ∈ PasXfis RPI set underK, then,KZf⊂V, see (17b), i.e.,

∀z_i|k ∈Zf, Kz_i|k ∈V, i≥N. (20) This indicates thatZf is RPI for any nominal system at any k≥0 under the controllerK, i.e.,

∀z_i|k ∈Zf, zi+1|k ∈Zf, i≥N. (21) Finally, we present a formal statement for the stability of the proposed LPVMPC and the recursive feasibility of its optimization problem (18).

Theorem 7 (Recursive Feasibility and Stability of the LPVMPC for Regulation): Consider the LPV system (1), suppose that Assumptions 4 and 5 are satisfied and let V_f in (18a) be a terminal cost satisfying the condition in Definition1 and the set Xf be a terminal set according to Definition2, then

(i) the optimization problem (18) is recursively feasible and (ii) the LPVMPC solution by (18) is asymptotically stabi-

lizing.

The proof of Theorem 7 follows the same lines as of the original approach in [17] taking into account the related simplifications considered here.

III. LPV MODELING OF A CLASS OF NL DYNAMICS In this section we present a systematic LPV modeling technique of a class of NL dynamics, which describes a wide range of mechanical systems. This is essential for the proposed MPC tracking problem in the next section.

(6)

Then, we apply the technique on the NL dynamics of a 2-DOF robotic manipulator as shown in SectionV-A.

Consider the Lagrangian formulation of the NL dynamics of ann_q-DOFs mechanical system given by

M(q)q(t)¨ +c(q,q)˙ +g(q)+τv(q)˙ =τ, (22) where q ∈ Rⁿ^q is the vector of generalized coordinates, M is the inertia matrix, which is assumed to be invertible, the vector cincludes Coriolis and centrifugal terms,g contains the terms derived from the potential energy, such as gravitational forces, τv denotes friction vector andτ is the vector of control inputs. The states of the system are commonly the vector

q^> q˙^>^>

. To simplify the derivation of a qLPV representation of NL systems given by (22), consider the the transformed states given as follows

x =

I 0 0 M(q)

q q˙

= q

M(q)q˙

, thus,

x˙ =

q˙ M(q)q¨+ ˙M(q)q˙

. (23)

Next, substituting the termM(q)q¨from (22) into (23), exploiting the relation between the terms M˙(q)q˙ andc(q,q), see,˙ e.g., [28], reformulating the vector g(q) in matrix form as shown below, and considering just viscous friction for the term τv(t), lead to the following continuous-time general qLPV state-space representation, which is equivalent to the NL model (22),

x˙ = ˜A(q,q)x˙ + ˜Bu, (24a)

y= ˜Cx (24b)

whereu=τ,

A(q˜ ,q)˙ =

0 A˜₁₂(q) A˜₂₁(q) A˜₂₂(q,q)˙

, B˜ =

0 Inq

, C˜ = In_q 0

(25) with A˜12(q) = M⁻¹(q),A˜21(q) is related to the termg(q) and A˜22(q,q) is related to the terms˙ c, τv. The scheduling parameter can be defined as a function of q and q, see˙ e.g., [29]. Moreover, the output of the system according to (24b) and (25) isq, which is often the controlled variable.

Note that only theAmatrix is parameter dependent, which is quite desired in several LPV control design approaches [6].

The structure of the system matrices in (25) will be very useful for developing the proposed MPC tracking problem in the next section.

Remark 8: In case of underactuated systems [30], which have fewer control inputs than DOFs, some elements of the vectorτ in (22) are dependent on each others or zeros, that might limit the controller authority on the system in comparison with that on fully actuated systems. Interestingly, the presented LPV modeling as in (24) can be used for underactuated systems; however, the elements dependency in τ should be

taken into account in the MPC optimization problem. That can be easily handled provided that such dependency is linear;

otherwise, nonlinear transformations can be performed onτ as commonly used in the context of nonlinear control [31].

Furthermore, since MPC is considered here to control the system, a discrete-time model should be obtained. For simplicity we use Euler’s forward (rectangular) discretization, which results in a discrete-time qLPV model as shown in (1) with

A(p)=I+ ˜A(q,q)T˙ _s, B= ˜BTs, C= ˜C, D=0 (26) whereA˜,B˜,C˜ are given in (25) andT_sdenotes the sampling time. The rectangular method is an approximative method of discretization; however, it has the important feature that it can preserve the linear dependence over the scheduling variables without introducing any extra complexity [2]. Furthermore, it preserves stability and yields a small discretization error provided that a suitable value ofTsis chosen, see [2] for more details.

IV. PRACTICAL LPVMPC FOR REFERENCE TRACKING Based on the LPV modeling presented in the previous section, we propose in this section a novel MPC formulation for LPV systems to track a given desired reference trajectory using the notion of admissible reference. For ensuring stability and recursive feasibility, a terminal cost as a tracking error penalty term is added to the stage cost of the related MPC optimization problem and a terminal constraint based on the concept of invariant set for tracking [23] is included; both are computed offline. The MPC optimization problem is formulated based on the nominal LTI model corresponding to the current value of the scheduling variable. Therefore, the control law is obtained online by solving a single QP problem of a complexity similar to that of LMPC. To deal with the uncertainty of the scheduling variable affecting the evolution of the state over the prediction horizon, we utilize the notion of rigid tubes presented in SectionIIusing a set Ssatisfying Assumption5. It can be computed offline taking into account the rate of variation of the scheduling parameter, which can reduce the size ofSleading to a low conservative design.

A. ADMISSIBLE REFERENCE

For the LPV system (1), any admissible steady state and input should satisfy the following equation

A(p)¯ −I B(p)¯ x¯

u¯

=0, p¯∈P (27) wherex¯andu¯are the steady state and input, respectively, and p¯is the steady-state value of the scheduling parameter, which represents —in case of qLPV models— the frozen scheduling parameter associated with (x¯,u). Note that (¯ x¯,u) belongs to¯ the null space of the left matrix in (27); moreover, provided that the system is controllable for allp¯ ∈P, the dimension of the null space isn_u.

(7)

Now, consider qLPV models of the nonlinear dynamics (22) as illustrated in SectionIIIand make use of the special structure in (26) with (25), we can rewrite (27) as

0 A˜₁₂(p) 0¯ A˜₂₁(p)¯ A˜₂₂(p)¯ I_n_q

x¯ u¯

=0, p¯ ∈P, (28) note here thatn_q=n_u, see Remark8. It holds that

x¯ u¯

=L(p)¯ ¯y_s, p¯∈P, (29) wherey¯_s ∈ Rⁿ^q characterizes the solution of (28) and L ∈ R⁽ⁿ^x⁺ⁿ^q^)×n^q can be given as

L(p)¯ =





 I_n_q

0

− ˜A21(p)¯





, (30) which spans the nullspace of the left matrix in (28). The parameter y¯_s represents the associated admissible output, which is related to x¯ by y¯_s = Cx, see (26) and (25). In¯ other words,y¯_s represents admissible references which can be tracked without steady-state tracking error. Moreover,p¯is related toy¯_s, this means that for anyy¯_sthere existsp¯function ofx¯for the qLPV models considered in SectionIII.

B. ADMISSIBLE RIST

Consider system (1) with the constraint sets X andUand the special matrices structure in (26) with (25), let also Assumption4be satisfied. Moreover, suppose that the system is controllable for allp∈Pby a control law as in (4), which meets the conditions in Definition 3. Therefore, the system matrixA(p)+BK is Schur for allp∈Pand the closed-loop system converges to (x¯,u) without violating the state and¯ input constraints provided that the initial state and (x¯,u)¯ belong to the corresponding robust invariant set for tracking.

This set can be considered as an RPI set for the augmented system

x(k+1) y¯_s

=

A(p(k))+BK −B(K1+ ˜A21(p))¯

0 In_q

| {z }

A¯^c(p(k),p)¯

x(k) y¯_s

(31)

subject to the system state and input constraints. The system (31) is derived by substituting (4) into (1a) taking into account (26), (25), using (29) and partitioningK =

K1 K2

such that K₁∈Rⁿ^u^×(n^x⁻ⁿ^u⁾. Note that for a given admissible¯y_sand the correspondingp, the state of the closed-loop system (31) will¯ converge to the relatedx¯according to (29) starting from any initial state and the associatedp ∈ P. Since the system (31) subjects to constraints, we want to construct an RPI set in which any initial state will converge to a steady state according toy¯s for anyp ∈ P. Then, that RPI set can be used as a terminal set for the proposed LPVMPC.

The set of constraints on (31) can be posed as follows X¯_λ= {(x,y¯s)| x∈X, ¯ys∈λY,¯ p¯∈P,

(Kx−(K1+ ˜A21(p))¯ y¯_s)∈U}, (32)

where ¯

Y=Projy(X) denotes the constraint set of the admissible output andλ is a scalar, for the moment, letλ = 1.

Now, define a set⊆ ¯X1as an admissible robust invariant set for tracking for the system (1) with (26) or as a robust positive invariant set for the augmented system (31), under the constraint set ¯

X1, if

∀ x

y¯s

∈: ¯A^c(p,p)¯ x

y¯s

∈, ⊆ ¯X1, ∀p∈P. (33) In order to attain the largest possible domain of attraction of the proposed MPC scheme, we should consider the maximum admissible RIST, which can be determined —using linear programming— similarly as computing a maximum RPI set for the system (31) constrained in ¯

X1. However, due to the unity eigenvalues of the augmented matrixA¯^cas shown in (31), the computation of maximum RPI set may not be finitely determined. Apparently, this can be resolved by considering the set ¯

Xλin (32) withλarbitrarily close to 1, note thatλ ∈ (0,1), see [32] for more details. This results in an admissible RIST set smaller but close to the maximum one.

Denote such a set as ¯

Xf ⊆ ¯Xλ, which will be employed for constructing the terminal set in the proposed MPC setting as shown below, therefore, ¯

Xfis the maximum⊆ ¯Xλ. C. OPTIMIZATION PROBLEM

Now, the proposed MPC optimization problem for reference tracking can be formulated as follows:

min v_0|k,· · ·,v_N−1|k y¯0|k,· · ·,¯y_N|k

N−1

X

i=0

kzi|k− ¯zi|kk²_Q+ kvi|k− ¯vi|kk²_R

+ k¯y_i|k−r_i|kk²

T + k¯y_N|k−r_N|kk²

T

+ ¯Vf(zN|k,¯zN|k) (34a)

subject to (16a-d), (34b)

z_N|k y¯_N|k

∈ ¯Zf (34c)

and the nominal system dynamics atp(k) withz_0|k =x(k), whereT = T^> 0,Q,Rare tuning matrices to achieve a desired tracking performance, for alli=0,1,· · · ,N,¯z_i|k,

¯v_i|kare nominal steady state and input, respectively, according to the desired reference trajectoryr_i|k, which are parameterized as in (29) by the corresponding nominal admissible outputy¯_i|k withp¯ = p(k),V¯f(z_N|k,¯z_N|k) is the terminal cost for tracking which is given by

V¯_f(z_N|k,z¯_N|k)= kz_N|k − ¯z_N|kk²

P, (35)

and satisfies the condition in Definition1and¯

Zfis a tightened terminal set given by

Z¯f= ¯Xf S×Proj_y(S)

, (36)

with ¯

Xfis an admissible RIST as discussed in SectionIV-B.

The nominal statesz_i|k are computed according to the nominal model in (19) given p0|k = p(k) (for computing the corresponding system matrix), where z0|k = x(k).

(8)

The tightened state and input constraint sets Z and V are computed offline from (17a) and (17b), respectively. The reference stepsr0|k,· · · ,vN|kare given in advance. The nominal control input movesv_0|k,· · ·,v_N−1|k, and the admissible output movesy¯_0|k,· · ·,y¯_N|kare the decision variables of the optimization problem. Finally, the values of¯z_i|k and¯v_i|k are substituted in terms ofy¯_i|k using (29) withp¯ =p(k) as these nominal steady state and input are related to the nominal model atp(k), i.e., (19), where its matrices computed atp(k) are frozen over the prediction horizon. Another way, one can use the nominal LPV model in (10) instead of (19); however, the system matrices as well as the matrix L in (30) should be updated at each step of the prediction horizon, which demands a slight increase in computations.

The main result of the paper is presented in the following theorem.

Theorem 9 (Recursive Feasibility and Stability of the LPVMPC for Tracking):Consider the LPV system (1) with the system matrices (26) and (25) and suppose that Assump- tions4and5are satisfied. Letrbe a given admissible reference trajectory,V¯_fin (35) be a terminal cost for the optimization problem (34) satisfying the condition in Definition1and Z¯fin (36) be its terminal set, where ¯

Xfsatisfies the invariance condition in (33) withand ¯

X1are replaced by ¯ Xfand ¯

Xλ, respectively, for a givenλ∈(0,1), then

(i) the optimization problem (34) is recursively feasible and (ii) the LPVMPC solution by (34) is asymptotically stabi-

lizing.

The proof of Theorem9is detailed in AppendixA.

Remark 10: The approach in [23] and [24] has considered fixed values for¯zi|k andv¯i|k for alliover N, thus, it handles piecewise constant references, moreover, that can lead to tracking with a large rise time. In contrast, the proposed approach is not restricted to that as the variables¯z_i|k andv¯_i|k are allowed to be changed overN via the admissible output moves y¯_0|k,· · ·,y¯_N|k, which does increase the number of decision variables, according to the given desired reference.

In this perspective, one might realizez¯_i|k andv¯_i|k as points on state and input trajectories related to the trajectory of the admissible reference rather than steady state and input values discussed in the previous section.

Remark 11: The termsk¯yi|k −ri|kk²

T,i= 0,1,· · · ,N in the cost function (34a) penalize the deviation between the nominal admissible output and the desired reference trajectory. This guarantees offset-free tracking provided that there exists an admissible output trajectory equal to the desired reference; otherwise, it steers the system to the closest admissible output.

Algorithm1summarizes the online implementation of the proposed approach. It can be executed provided that Assump- tions 4(ii) is satisfied and given the reference trajectory r, the matricesP,Q,R,T and the setsX,Z,U,Vand ¯

Zfwhich can be computed offline as shown in the next section. At every samplek, the value ofp(k) is measured and used to update the matrixA using (26), (25), and hence the block matrixA˜21,

and the matrixLas in (30). Then, givenx(k) the optimization problem (34) is solved and the receding horizon concept is used, where the control sampleu(k) =v0|k is applied to the system.

Algorithm 1The Proposed LPVMPC for Tracking Require: r,P,Q,R,T,X,Z,U,Vand¯

Zf FOffline Initializationx(0),p(0) andk=0

Repeat FOnline

1: Measurep(k) and updateA(p_0|k) andL(p_0|k).

2: Measurex(k) and solve the optimization (34).

3: Implement the control sampleu(k)=v0|k.

4: k←k+1

D. OFFLINE COMPUTATIONS

The proposed LPVMPC approach involves offline computations including the cost functionV¯_f and the sets Z,V, ¯

Zf. ComputingV¯ftogether with the controller gainK according to Definition 1 is a standard LMI problem, see [33] for more details. Regarding the computations of the tightened constraint sets Z, V, ¯

Zf, it is based on the setsX, U, ¯ Xf, respectively, as well as the set S as shown in (17a,b) and (36). Determining the set ¯

Xfis also a standard problem using linear programming which can be computed as a maximum RPI set, as shown in [26], for the system (32) subjected to the constraint set ¯

X_λwithλchosen close to 1 so that ¯

Xfis finitely determined. Concerning the setSsatisfying Assumption5, we propose two ways as follows:

1) One way is to considerS a minimal robust positive invariant (mRPI) set [34] based on a disturbance set given asCo{(A(p)−A0)x+(B(p)−B0)u| p∈P,x∈ X,u∈U}corresponding to the whole scheduling range P, where (A₀,B₀) are the related nominal system matrices evaluated at the center ofP. In this way, we ignore the fact that the rate of variation ofpis bounded. The approach proposed in [34] can be employed to compute an outer approximation of the mRPI set.

2) Another way to computeS, letN and the bound on the rate of variation ofp, i.e., dp^max, be given, calculate the set S_N|k using (13) on a grid points of P, which we denote asS_N|kⁱ^g at every grid pointig =1,· · ·,ng, whereng is the number of the grid points. Next, use these sets to compute

S=α·Co







ng

[

ig=1

S_N|kⁱ^g







, (37)

whereα≥1. Note that, any of the setsS_N|kⁱ^g is a PC-set according to the formulation in Section II. Finally, the obtainedScan be verified on a denser grid to check the validity of Assumption5; otherwise, it should be enlarged usingαand verified again.

The first way is guaranteed to verify condition (15); however, it might be overly conservative, especially, when dp^max

(9)

is relatively small or N is short, resulting in a significant reduction in the size ofZftogether with the domain of attraction, then the second method could provide a better solution provided that condition (15) is verified on very dense grid points. It is recommend to compare the obtainedSfrom both approaches.

V. APPLICATION TO THE CRS A465 ROBOT

In the following, we present the results of implementing the proposed LPVMPC for tracking on the CRS A465 Robotic manipulator shown in Fig2. It has six rotational joints actuated by DC motors. The angular displacements of the motor shafts are measured by incremental encoders. In this work, we consider the shoulder and the elbow corresponding toq1

andq₂, respectively, they are the most challenging links to control, since they are affected by gravity, inertial, centripetal, Coriolis and friction torques; the other links are fixed during the experiments. To demonstrate the quality of the proposed control approach we consider two practical trajectories [35]

to be tracked and we compare with other MPC scheme developed recently in [22].

FIGURE 2. The CRS A465 robotic manipulator and a side view of the 2-DOF model.

A. LPV MODELING

Based on the formulation introduced in SectionIII, we derive now a qLPV model for a 2-DOF of the CRS A465 robotic manipulator. According to (22),q,q,˙ q¨ ∈ R² are the joint angular positions, velocities and acceleration, respectively, and the nonlinear matrices are given by

M =

b7+b2+2b3cos(q2) b2+b3cos(q2) b7−b8+b3cos(q2) b7

, c=

−b3q˙²₂sin(q2)−2b3q˙1q˙2sin(q2) b₃q˙²₁sin(q₂)

, g =

−b4sin(q1+q2)−b5sin(q1)

−b4sin(q1+q2)

, τv= b6q˙1

b9q˙2

(38) where b1,· · · ,b9 are grouped parameters related to the dynamic and kinematic parameters of the robot. We use the parameters of [4], which were experimentally identified for the same robot considered here, the parameters are shown in Table 1. The physical limits of the robot movements are specified by the manufacturer of the robot as the following bounds on its angular positions and velocities

|q1|≤90^◦, |q2|≤110^◦, | ˙q1|≤180^◦, | ˙q2|≤180^◦, (39)

TABLE 1.Constant parameters in(38).

and its input constraints are

|τ1|≤5, |τ2|≤5, (40) which are based on the actuators limits.

Therefore, the qLPV representation of the system is given by (24) with (25), where

A˜12(q)=M⁻¹(q) (41a)

A˜21(q)= b₄sin(q₁+q₂) q1+q2

I2+b₅sin(q₁) q1

1 0

0 0

(41b) A˜₂₂(q,q)˙ = −

b₆ 0 0 b₉

+b₃q˙₁sin(q₂) 0 0

1 1

M⁻¹(q). (41c) To further simplify the obtained qLPV model, we neglect the second-order termb3q˙1sin(q2) in (41c), this gives

A˜₂₂(q)= −

b₆ 0 0 b₉

M⁻¹(q), (42) therefore, the system matrixA˜ becomes dependent only on q, and hence, measuring A˜ in practice is independent of derivatives of signals. In general, the error due to such sim- plification is usually small, unless when bothq₂andq˙₁are very close to their limits, which is an irregular scenario. Next, we define the scheduling variables of the qLPV model based on (41a,b) and (42), which indicate 4 variables yieldingA˜ affine dependent matrix; they are shown together with their bounds and the bounds on their derivatives according to (39) in Table2. Finally, the discrete-time model of the system can be determined as in (26) taking into account the blocks ofA˜ as shown in (41a,b) and (42).

TABLE 2.Scheduling parameters data.

B. SIMULATION RESULTS

We consider the discrete-time qLPV model developed in SectionV-A, with sampling timeT_s =0.01s, which compro- mises the execution time of the MPC algorithm as well as the bandwidth of the system. Moreover, we scale the scheduling variables in Table2to be within±1, which is preferred for numerical reasons. Note that scalingpchanges the constant matrices associated with the affine dependency onp, i.e., the

(10)

matricesA^jin (5). In discrete-time, the rate of change ofpis defined as dp= ˙pTsand its bounds are computed accordingly using (39) and the definition ofpin Table2, then, they are redefined based on the scaledpas follows:

|dp₁| ≤0.073, |dp₂|≤0.117

|dp₃| ≤0.086, |dp₄|≤0.078. (43) To compute the terminal costV¯_f(35), the matrixP∈R^4×4 together with the associated controllerK ∈R^2×4are obtained using LMIs;K is required for calculating the sets ¯

XfandS. The RIST ¯

Xfhas been computed as discussed in SectionIV withλ = 0.975, which can lead to reasonably offset free tracking ranges of

|q₁|≤87.75^◦, |q₂|≤107.25^◦, (44) c.f., (39); note that larger value than thatλcould not allow convergence while computing the set ¯

Xf. To calculate the set Ssatisfying Assumption5, a prediction horizonN =4 has been chosen and used afterward in the online implementation, then, the second method in Section IV-D has been carried out to computeS, which is a reasonable choice here due to the small values of dp_max as shown in (43) and the short prediction horizon. According to (37), we obtained the setS withng=3⁴andα=1.1 and it has been verified on a very dense grid of 81⁴points on the set P. Next, we performed set tightening to obtain the constraint setsZ,VandZ¯fusing (17a), (17b) and (36), respectively, which are the required sets to implement Algorithm1. It is worth mentioning that all the involved set computations and operations have been performed using the Multi-Parametric Toolbox 3.0 [36]. To solve the optimization problem (34), the weighting matrices for the stage cost have been tuned toQ = 10·(1,1,0,0), R=0.001·diag(1,1) andT =1000·diag(1,1).

To evaluate Algorithm1, it has been implemented in simulation to track a sinusoidal reference for bothq₁andq₂which covers the whole range of both as given in (39) and results in an ellipse like shape reference trajectory in the Cartesian space. A projection of the state-phase plane on the (x₁-x₂)- plane is shown in Fig. 3, note thatx1 = q1 andx2 = q2

according to (23). The figure demonstrates the convergence of the state trajectories ofq1andq2to the admissible command trajectory based on the smaller range (44), c.f., (39).

This leads to a small steady-state tracking error with respect to the desired Cartesian command trajectory. In Fig.3 it is also shown the projection of the computed set S onto the (x₁, x2)-plane, which is quite very small leading to a very small difference between the state constraint set Xand the tightened one Z as depicted in the figure as well by their projections onto the (x₁,x₂)-plane. In addition, the projection of the terminal set¯

Zfonto the (x₁,x₂)-plane is shown in Fig.3.

Moreover, for comparison, we have computed the terminal set Z⁰_f corresponding to the steady state at the origin, its projection onto the (x₁,x₂)-plane is depicted in Fig.3. Obviously, it is much smaller than the terminal set based on the proposed RIST, which implies larger domain of attraction. For further

FIGURE 3. Phase plane on the (x₁,x₂)-plane (the transformed state as shown in (23)): the reference (yellow) and the simulated state trajectories of (x₁,x₂) (dashed-blue), and the projection of the setsX(sold-gray), Z(dashed-black),Z¯_f(filled-light-gray),Z⁰_f (filled-gray) andS(the smallest set) onto the (x₁,x₂)-plane. Note thatx₁=q₁andx₂=q₂.

FIGURE 4. Phase plane on the (x₃,x₄)-plane (the transformed state (23)):

the projection of the setsX(sold-gray),Z(dashed-black),Z¯f (filled-light-gray),Z⁰_f (filled-gray) andS(the smallest set) onto the (x₃,x₄)-plane.

illustration, Fig.4shows the projections of the setsX,Z, ¯ Zf, Z⁰_f andSonto the (x₃,x₄)-plane.

C. EXPERIMENTAL RESULTS

To evaluate the proposed LPVMPC approach in real-time, we consider here two tracking scenarios: In the first one, the robot tracks typical trajectories from practice, whereas, in the second scenario the robot tracks command trajectories beyond its work space. For comparison we consider the MPC approach in [22]. It is an NMPC scheme based on qLPV models for NL systems with stability guarantees, where the scheduling parameters are predicted in an iterative procedure over the prediction horizon without convergence guarantees and the online nonlinear optimization problem is solved as a sequence of QPs, which might be computationally demanding in comparison with the proposed approach. Moreover,