# Dynamic Output Feedback Robust Model Predictive Control via Zonotopic Set-Membership Estimation for Constrained Quasi-LPV Systems

Authors: Xubin Ping; Ning Sun
Journal: Journal of Applied Mathematics
 (2015)

---

## Abstract

For the quasi-linear parameter varying (quasi-LPV) system with bounded disturbance, a synthesis approach of dynamic output feedback robust model predictive control (OFRMPC) is investigated. The estimation error set is represented by a zonotope and refreshed by the zonotopic set-membership estimation method. By properly refreshing the estimation error set online, the bounds of true state at the next sampling time can be obtained. Furthermore, the feasibility of the main optimization problem at the next sampling time can be determined at the current time. A numerical example is given to illustrate the effectiveness of the approach.


---

## Body

## 1. Introduction

Model predictive control (MPC) or receding horizon control is a class of optimization based control methods, which explicitly utilizes process model parameters and measurements to optimize a performance function. At the current time, the optimization problem is carried out and then generates the optimal control input sequence. However, only the first element of control inputs is implemented. At the next sampling time, the output measurement is updated and the prediction horizon is shifted one step forward, and then the optimization problem is repeated. Due to its ability to handle the systems with hard constraints, it has attracted many researchers’ attention, from both academic community and industrial society, for example, [1–4]. However, in real processes, the precise model parameters are seldomly available. Hence, robust MPC (RMPC) is more practical for real applications.Linear parameter varying (LPV) systems are the systems whose parameters take their values in the prespecified sets and the dynamic characteristics depend on the time-varying parameters. When the time-varying parameters of LPV systems can be exactly known at the current time, the systems are called quasi-LPV (quasi-linear parameter varying) systems. It means that the model parameters are exactly known at the current time, but their future evolutions are uncertain and contained in the prescribed bounded sets. In RMPC, the model parametric uncertainty can be dealt with within the frame of LPV systems. For the online RMPC, at each time, a min-max optimization problem is often utilized to minimize the performance function of LPV systems, which considers all the possible realization of model parametric uncertainty [5]. Reference [6] extends the approaches in [5] to the quasi-LPV system with a quasi-min-max optimization. In [7, 8], the controller designment considers the time-varying parameters of quasi-LPV systems having the bounds on their rate of variation. Reference [9] considers the quasi-LPV systems with a parameter-dependent control law.However, in [5–9], the true states are assumed to be known and the bounded disturbance is not considered. In real processes, the true states are usually unmeasurable, and only the outputs with disturbance are available. In these situations, output feedback RMPC (OFRMPC) is more practical than the state feedback one for real applications. For OFRMPC with bounded disturbance and without model parametric uncertainty, see [10, 11]. For OFRMPC with model parametric uncertainty and without bounded disturbance, see [12–14]. For OFRMPC with both polytopic uncertainty and bounded disturbance, see [15–20].For OFRMPC, the bounds of estimation error set represent a kind of uncertainty, which has to be considered in the robust stability and physical constraints. When the model parametric uncertainty, bounds of estimation error set, and disturbance are known and contained in the bounded sets, the evolution of state estimation is a compact set which is consistent with the bounded uncertainties. The set-membership estimation considers nonstatistical description of model parameters, estimation error sets, and disturbance in the form of bounded uncertainties. When the set-membership estimation is utilized to state estimation, an appropriate representation of estimation error sets should consider the tradeoff between the precise representation of estimation error sets and the online computational burden. Polytopes can be used to represent exactly the estimation error set. However, if the vertices of polytopes increase dramatically, the computational burden becomes quickly prohibitive [21]. Ellipsoids can overcome the drawback by efficient computation. However, the Minkowski sum of ellipsoids is not necessary an ellipsoid, and an outer approximation ellipsoid is often used to estimate the Minkowski sum of ellipsoids, which may lead to the conservatism of state estimation. In recent years, zonotopes (a particular class of polytopes) have received more attention because of their advantages of having better precision in comparison with ellipsoidal sets and less complexity compared to polyhedral sets. Hence, several set-membership estimation methods based on zonotopes [22–25] have been proposed.The present paper considers a synthesis approach of dynamic OFRMPC for the quasi-LPV system with bounded disturbance. The main optimization problem is similar to that in [16, 19], while the estimation error sets are represented by zonotopes. The main contribution is a combination between the main optimization problem that calculates control parameters and the set-membership estimation based on zonotopic computation. By properly refreshing the estimation error sets, it can obtain the precise bounds of true states at the next sampling time. Furthermore, the feasibility of the main optimization problem at the next sampling time can be determined by checking a feasible problem at the current time.


## 2. Notations, Basic Definitions, and Properties

For any vectorx and positive-definite matrixW,xW2≜xTWx.x(i∣k) is the value of x at time k+i, predicted at time k. I is the identity matrix with appropriate dimension. εM≜ξ∣ξTMξ≤1 denotes the ellipsoid associated with the symmetric positive-definite matrix M. All vector inequalities are interpreted in an element-wise sense. An element belonging to Co· means that it is a convex combination of the elements in ·, with the scalar combing coefficients nonnegative and their sum equal to 1. The symbol “⋆” induces a symmetric structure in the matrix inequalities. A value with superscript “*” means that it is the optimal solution of the optimization problem. The time-dependence(k)of the MPC decision variables is often omitted for brevity.Let an interval[a;b] be defined as the set {x:a≤x≤b}. The unitary interval is B=[-1;1]. A box ([a1;b1],…,[an;bn])T is an interval vector. A unitary box in Rm, denoted by Bm, is a box composed by m unitary intervals. The Minkowski sum of two sets X and Y is defined by X⊕Y={x+y:x∈X,y∈Y}.Definition 1.
Am-zonotope in Rn can be defined as the linear image of a m-dimensional hypercube in Rn, where m is the order of the zonotope and satisfies m≥n. Given a vector p∈Rn and a matrix H∈Rn×m, a m-zonotope is the following set:(1)Z=p⊕HBm=x∈Rn∣x=p+Hz,z∈Bm.This is the Minkowski sum of them-segments defined by m columns of matrix H in Rn. The center of zonotope Z is vector p; H=h1,…,hi,…,hm, with hi∈Rn,  i∈{1,…,m}, are the generators of the zonotope. Some properties for zonotopic computation are given as follows.Property 1 (the Minkowski sum of two zonotopes [22]).
Given two zonotopesZ1=p1⊕H1B1m∈Rn and Z2=p2⊕H2B2m∈Rn, the Minkowski sum of two zonotopes is also a zonotope defined by Z=Z1⊕Z2=(p1+p2)⊕[H1H2]Bm1+m2∈Rn.Property 2 (the image of a zonotope [22]).
The image of a zonotopeZ=p⊕HBm∈Rn by a linear mapping K can be computed by a standard matrix product KZ=(Kp)⊕(KH)Bm∈Rn.Property 3 (zonotope reduction [22, 23]).
Given a zonotopeZ=p⊕HBm∈Rn and integer f with n<f<m, denote by H^ the matrix resulting from the reordering of the columns of matrix H in decreasing order of Euclidean norm and by O∈Rn×n a diagonal matrix, satisfying the following conditions:(2)H^=h^1,…,h^i,…,h^m,withh^i2≥h^i+12,(3)Oii=∑j=f-n+1mH^ij,withi=1,…,n,and then Z⊆Z^≜p⊕[H^TO]Bf∈Rn, where H^T is obtained from the first f-n columns of matrix H^.Remark 2.
From Property3, it can obtain that Z^ is an outer approximation of Z with a lower order, which can be utilized to estimate a zonotope with high order. In Section 4.2, by applying Property 3, the order of the estimation error set can be reduced.Zonotopes are the special case of polytopes and can be represented by half-space representation (H-representation) and vertex representation (V-representation). The definitions are as follows.Definition 3 (H-representation of a polytope).
Forq half-spaces, a convex polytope P is the set: (4)P=x∈Rn∣Gx≤g,G∈Rq×n,g∈Rq×1.Definition 4 (V-representation of a polytope).
A polytopeP is the set constructed by the convex combination of finite vertices χi∈Rn,  i∈{1,…,VP}; that is,(5)P=Coχ1,…,χVP=∑i=1VPαiχi∣χi∈Rn,αi∈R,αi≥0,∑i=1VPαi=1,where VP is the total number of vertices related to the V-representation of polytope P. From Definition 1, the unitary box Bm in Rm has 2m vertices, and each vertex of Bm is a m-dimensional vector with the element being 1 or -1. Hence, according to (1), at most 2m vertices of zonotope Z can be calculated. We define the vertex which is the interior point of Z as redundant vertex. To represent the zonotope Z by V-representation through VZ vertices, a Quickhull Algorithm in [26] can be employed to remove the redundant vertices. When the dimension n of zonotope Z is fixed, the relationship between VZ and m can be obtained; for example, when n=2, VZ=2m; n=3, VZ=m(m-1)+2, respectively [27]. For high-dimensional systems, although the relationship between VZ and m is complex, the upper bound of VZ can be calculated in polynomial time for the fixed dimension of zonotopes [28].Figure1 shows a zonotope Z-=p⊕HB3, where p=65, H=2.5-20.50-1-2, m=3; the three generators are h1=2.50, h2=-2-1, h3=0.5-2, which are the three columns of matrix H, respectively [27]. Based on Definition 1, the vertices of zonotope Z- can be calculated for the different 23 vertices of unitary box B3, where the different combinations of “+” and “−” denote the corresponding vertices of Z-. In order to describe the zonotope Z- by V-representation through 6 (Vz=2×3) vertices, the redundant vertices denoted by “(-,-,+)” and “(+,+,-)” are removed.Figure 1
Representation of a zonotopeZ- by generators and vertices.


## 3. Problem Statement

Consider the following discrete-time uncertain LPV system:(6)xk+1=Akxk+Bkuk+Dkwk,yk=Ckxk+Ekwk,where u∈Rnu, x∈Rnx, y∈Rny, and w∈Rnw are input, state, output, and disturbance, respectively.We make the following assumptions on system (6).(A1)
The time-varying system matrices[A∣B∣C∣D∣E](k) are known to vary within a polytope Ω, that is, [A∣B∣C∣D∣E](k)∈Ω≜Co{[Al∣Bl∣Cl∣Dl∣El], l={1,…,L}}, and there exist nonnegative coefficients λl(k), where λl(k)’s are exactly known at current time k, such that ∑l=1L‍λl(k)=1 and [A∣B∣C∣D∣E](k)=∑l=1L‍λl(k)[Al∣Bl∣Cl∣Dl∣El].(A2)
The bounded disturbance satisfiesw(k)∈εPw, and an outer approximation of the bounded ellipsoidal disturbance set can be represented by a zonotope D; that is,(7)εPw⊆D≜0⊕HwBd,Hw∈Rnw×d,d≥nw,where matrix Hw is constructed by d generators of zonotope D.(A3)
The input and output constraints are required to satisfy(8)-u-≤uk≤u-,-ψ-≤Ψyk+1≤ψ-,where u-=[u-1,u-2,…,u-nu]T; u-j>0, j∈{1,…,nu}; ψ-=[ψ-1,ψ-2,…,ψ-q]T; ψ-j>0, j∈{1,…,q}; Ψ∈Rq×ny.For the system (6), the dynamic output feedback controller is of the following form [16, 19]: (9)xci+1∣k=Aci∣kxci∣k+Bci∣kyi∣k,ui∣k=Cci∣kxci∣k+Dci∣kyi∣k,where xc∈Rnx, {Ac(i∣k),Bc(i∣k)} are the controller gain matrices, and {Cc(i∣k),Dc(i∣k)} are the feedback gain matrices. The controller parameters {Ac(i∣k),Bc(i∣k),Cc(i∣k),Dc(i∣k)} take the parameter-dependent form as(10)Aci∣k=∑l=1L∑j=1Lλlk+iλjk+iA-clj,Dci∣k=D-c,Bci∣k=∑l=1Lλlk+iB-cl,Cci∣k=∑j=1Lλjk+iC-cj.Hence, the augmented closed-loop system, based on (9) and the predictions made by (6), is (11)x~i+1∣k=Φi,kx~i∣k+Γi,kwk+i,∀i≥0,x~0∣k=x~k,where (12)x~=xT,xcTT,Φlj=Al+BlD-cCjBlC-cjB-clCjA-clj,Γlj=BlD-cEj+DlB-clEj,Φi,k∣Γi,k=∑l=1L∑j=1Lλlk+iλjk+iΦlj∣Γlj.At each timek, the true state x(k) is unmeasurable. Hence it is necessary to use its bounds to represent the true state. Suppose, at time k, x(k) is confined by a m(k) order zonotope represented by(13)xk∈Xk=xck⊕HkBmk,kkkHk∈Rnx×mk,mk≥nx,and then by applying Definition 1 and removing the redundant vertices, the true state x(k) can be represented by the following convex set:(14)xk∈Coχ1k,χ2k,…,χVXkk,where χr(k),  r∈{1,2,…,VX(k)}, are the vertices of X(k) and  VX(k) is the total number of vertices related to the V-representation of zonotope X(k) at time k.


## 4. Synthesis Approach of Dynamic OFRMPC

### 4.1. The Main Optimization Problem

In the dynamic OFRMPC approach, at each timek, based on the dynamic OFRMPC approaches in [16, 19], the following main optimization problem is considered:(15)min⁡γ,Q,A-clj,B-cl,C-cj,D-cmax⁡A∣B∣C∣D∣Ek+i∈Ω,wk+i∈εPwγ,(16)s.t.x~k∈εQ-1,(17)x~i∣kQ-12≥1⟹x~i∣kQ-12-x~i+1∣kQ-12≥1γyi∣kQ2+ui∣kR2,∀i≥0,(18)-u-≤ui∣k≤u-,∀i≥0,(19)-ψ-≤Ψyi+1∣k≤ψ-,∀i≥0,where γ is the performance cost, (16) is the augmented state constraint, (17) is the stability/optimality condition, (18) is the input constraint, (19) is the output constraint, the matrices Q, R are the positive-definite weighting matrices, and Q satisfies the following conditions:(20)Q-1=M=M1M2TM2M3,Q=Q1Q2TQ2Q3.Choosing M2=-M1, then M3=M1(M1-Q1-1)M1, Q2=Q3=Q1-M1-1.In the following, some previous results mainly in [16, 19] are reviewed to deal with constraint conditions (16)–(19). By applying the definition of augmented state and (14), constraint (16) holds if there exist positive-definite matrices {M1,Q1,M-3,M4} such that(21)χrkxckTM1-M1-M1M-3χrkxck≤1,kkkkkkkkkkkkikr∈1,2,…,VXk,(22)M-3-M1-M1M4≥0,M1-M4IIQ1≥0.For all x~(i∣k) and all admissible w(k+i), condition (17) holds if there exist a scalar α and matrices {A^clj,B^cl,C^cj,D^c} such that(23)ΥQBi,k=∑l=1L∑j=1Lλlk+iλjk+iΥljQB≥0,∀i≥0,ΥljQB=1-αM1⋆⋆⋆⋆⋆⋆1-αI1-αQ1⋆⋆⋆⋆⋆00αPw⋆⋆⋆⋆Al+BlD^cCjAlQ1+BlC^cjBlD^cEj+DlQ1⋆⋆⋆M1Al+B^clCjA^cljB^clEj+M1DlIM1⋆⋆Q1/2CjQ1/2CjQ1Q1/2Ej00γI⋆R1/2D^cCjR1/2C^cjR1/2D^cEj000γI,with (10) being parameterized as(24)D-c=D^c,C-cj=D-cCjQ1-C^cjM1-1-Q1-1,B-cl=M1-1M1BlD-c-B^cl,A-clj=M1-1A^clj-AlQ1-BlD-cCjQ1+B-clCjQ1×M1-1-Q1-1+BlC-cj.Suppose, at timek, constrains (21), (22), and (23) are satisfied; there exist positive-definite matrices {U,Ξ~} such that(25)Υju=M1⋆⋆⋆IQ1⋆⋆00Pw⋆2D^cCj2C^cj2D^cEjU≥0,kikUss≤u-s2,j=1,…,L,s=1,…,nu,Υhi,k=∑l=1L∑j=1Lλlk+iλjk+iΥljh≥0,kkkkkkikh=1,…,L,∀i≥0,Υljh=M1⋆⋆⋆IQ1⋆⋆00Pw⋆ΨChAl+BlD^cCjΨChAlQ1+BlC^cjΨChBlD^cEj+DlΞ~,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkΞ~ss≤ψ~s2,s=1,…,q,ψ~s=12ψ-s-ζs,ζs=max⁡hΨEhPw-1EhTΨTss,where Uss(Ξ~ss) is the sth diagonal element of U(Ξ~), and then by applying (10) and (24), the constraints (18) and (19) are satisfied.From the above conditions, it is shown that the main optimization problem (15)–(19) can be approximated by(26)min⁡γ,α,Q1,M1,M-3,M4,A^clj,B^cl,C^cj,D^c,U,Ξ~γ,is.t.21,22,23,25.At each time k, if the main optimization problem (26) is solved, then (11) is quadratically bounded [29, 30] within a common Lyapunov matrix Q-1, and the controller parameters are calculated by (10) and (24). The main optimization problem (26) can be solved by an LMI (linear matrix inequality) tool, which is the same as in [16, 19].Remark 5.
In [16], the true state constraint in the main optimization problem is represented by the V-representation of general polyhedral sets. In [19], to guarantee the recursive feasibility of the optimization problem based on polyhedral bounds, the true state constraint in the main optimization problem utilizes either the ellipsoidal set or the polyhedral set, with the latter being used if and only if it is contained in the former. Compared with [16, 19], the present paper utilizes the V-representation of zonotopes to represent the bounds of true state.


### 4.2. The Refreshment of Estimation Error Sets

Define the estimation errore(k)≜x(k)-xc(k). At each time k, for the quasi-LPV system (6), the model parameters [A∣B∣C∣D∣E](k)are exactly known at the current time k. If the controller parameters {Ac(k),Bc(k),Cc(k),Dc(k)} are obtained, then according to (11) and the definition of estimation error, it is shown that(27)ek+1=∑l=1L∑j=1Lλlkλjk×Al+BlD-cCj-B-clCjxk+BlC-cj-A-cljxck+BlD-cEj+Dl-B-clEjwk=Ak+BkDckCk-BckCkxk+BkCck-Ackxck+BkDckEk+Dk-BckEkwk.By applying (7) and (13) and Properties 1 and 2, it can be obtained that(28)ek+1∈Ek+1≜Ak+BkDckCk-BckCk+BkCck-Ackxck⊕Ak+BkDckCk-BckCkHk,BkDckEk+Dk-BckEkHw×Bmk+d.From (28), the estimation error set at time k+1 can be represented by a [m(k)+d]-zonotope E(k+1)=p(k+1)⊕H-(k+1)Bm(k)+d, where(29)pk+1=Ak+BkDckCk-BckCk+BkCck-Ackxck,H-k+1=Ak+BkDckCk-BckCkHk,BkDckEk+Dk-BckEkHw.It can easily be obtained from (28) that, with the evolution of time, the order of zonotope E(k+1) will increase. Accordingly, the generators and the vertices of zonotope E(k+1) will increase, which may lead to more LMIs in the main optimization problem (26) to enforce the augmented state constraints and then increase the online computational burden of the main optimization problem (26). To decrease the online computational burden of the main optimization problem (26), we impose the restriction on the order of E(k+1) by m(k)+d≤Nz. When m(k)+d>Nz, Property 3 is applied to construct a low order Nz-zonotope E^(k+1)=p(k+1)⊕H^(k+1)BNz for outer approximation of the high order zonotope E(k+1); that is, E(k+1)⊆E^(k+1), where(30)H^k+1=H^Tk+1Ok+1,kkikkH^k+1∈Rnx×Nz,and H^T(k+1)∈Rnx×[Nz-nx] and O(k+1)∈Rnx×nx are calculated according to (2) and (3), respectively. From the above derivations, it can be obtained that the estimation error set at time k+1 is one of the following forms:(31)Pk+1=Ek+1,Hk+1=H-k+1,mk+1=mk+d,whenmk+d≤Nz;Pk+1=E^k+1,Hk+1=H^k+1,mk+1=Nz,whenmk+d>Nz.When the estimation error setP(k+1) is obtained, according to Definition 1, at most 2m(k+1) vertices can be calculated. By removing the redundant vertices, the set P(k+1) can be described by V-representation through VP(k+1) vertices. The vertices of the true state set X(k+1) are the same as P(k+1); that is, VX(k+1)=VP(k+1). Hence the V-representation of the true state at time k+1 is(32)xk+1∈Coχ1k+1,χ2k+1,…,χVXk+1k+1,χrk+1=xck+1+vrk+1,kkir∈1,2,…,VXk+1.Remark 6.
In Section2, we have mentioned the relationship between Vz (the number of vertices related to zonotopes), m (the order of zonotopes), and the dimension of zonotopes. For the refreshment of estimation error sets in Section 4.2, when the dimension of true state is nx=2, VX(k+1)=VP(k+1)=2m(k+1); when nx=3, VX(k+1)=VP(k+1)=m(k+1)[m(k+1)-1]+2, respectively [27]. For high-dimensional systems, when nx is fixed, the upper bound of VX(k+1) (or VP(k+1)) is related to m(k+1) and can be calculated in polynomial time [28]. In [16], VX(k+1) (or VP(k+1)) only depends on the dimension of true state; that is, VX(k+1)=VP(k+1)=2nx. In [19], the upper bound of VX(k+1) (or VP(k+1)) is bounded by properly selecting an integer p^.
In [16, 19], the outer approximation of ellipsoidal bounded disturbance is represented by the V-representation of the polytope. The present paper utilizes zonotopes to outer approximation of εPw. The outer approximation of ellipsoids introduces conservatism. However, by properly increasing the vertices of outer approximation polyhedral sets in [16, 19] and choosing matrix Hw with proper generators, the conservatism can be reduced.


### 4.3. The Feasibility of the Optimization Problem at the Next Sampling Time

The original optimization problem (15)–(19) is recursive feasible; that is, it is feasible for all k>0 if and only if it is feasible at time k=0. However, by using the bounds of true state to handle the uncertainty of true state, the recursive feasibility may be lost for the main optimization problem (26). If the main optimization problem (26) is solved at time k, then the optimal solutions {α,γ,A^clj,B^cl,C^cj,D^c,Q1,M1,U,Ξ~}*(k) are obtained, and x~(k)∈ε[Q*(k)]-1. At time k+1, by choosing(33)M4k+1=M1*k-Q1*k-1,M1k+1=M1*k,Q1k+1=Q1*k,M-3k+1=M3k+1=M1*kM1*k-Q1*k-1-1M1*kT,then (22) is satisfied as follows:(34)M-3k+1-M1k+1-M1k+1M4k+1≥0,M1k+1-M4k+1IIQ1k+1≥0.Based on (33) and (34) and the definitions of Q and M, then, by applying the Schur complement, it can be obtained that if the problem (35) is feasible:(35)findη,s.t.η⋆⋆χrk+1Q1*k⋆xck+1Q2*kQ3*k≥0,mmi0<η≤1,r∈1,2…,VXk+1,then (21) is satisfied and x~(k+1)∈ε[Q*(k)]-1. In the main optimization problem (26), only (21) and (22) are related to the process state information. By further choosing {α,γ,A^clj,B^cl,C^cj,D^c,U,Ξ~}(k+1)={α,γ,A^clj,B^cl,C^cj,D^c,U,Ξ~}*(k), the feasible solutions for the main optimization problem (26) at timek+1can be constructed. The feasible solutions for (21) at timek+1can be directly constructed through (33). Hence, the feasibility of (35) determines not only whether x~(k+1) is contained in the ellipsoidal set ε[Q*(k)]-1 or not, but also the feasibility of the main optimization problem (26) at time k+1.Remark 7.
In [16], the refreshment of estimation error sets and the feasibility of the main optimization problem at the next sampling time are performed simultaneously by an iterative auxiliary optimization problem. Furthermore, it involves the H-representation of estimation error set in the auxiliary optimization problem. In [19], to guarantee the recursive feasibility of the main optimization problem based on polyhedral sets, the recursions of both polyhedral sets and ellipsoidal sets are obtained. At the next sampling time, the polyhedral set is used in the main optimization problem if and only if it is contained in the ellipsoidal set. The present paper refreshes the estimation error set via zonotopic computations and is represented by V-representation. Then problem (35) is applied to determine whether the main optimization problem (26) is solved or not at the next sampling time.


### 4.4. Overall Solutions

Algorithm 8.
Choose a feasiblexc(0) and the initial constrain of x(0) represented by a zonotope xc(0)⊕H(0)Bm(0), where H(0)∈Rnx×m(0). Let pro=1. Select the maximal order of zonotope related to the estimation error set as Nz. At each time k≥0, perform the following steps.(1)
Ifpro=1, then solve the main optimization problem (26) and obtain the controller parameters by (10) and (24).(2)
Ifpro=0, then inherit the previous controller parameters.(3)
Implementu(k)=Cc(k)xc(k)+Dc(k)y(k), and calculate xc(k+1)=Ac(k)xc(k)+Bc(k)y(k).(4)
CalculateE(k+1) by (28). If m(k)+d≤Nz, then let P(k+1)=E(k+1); else then calculate E^(k+1) by (30) and let P(k+1)=E^(k+1).(5)
Check the feasibility of problem (35). If problem (35) is feasible, then letpro=1; else then let pro=0.The complexity analysis of Algorithm8 and the related algorithms for worst cases are listed in Table 1. For the main optimization problem in Table 1, an LMI tool is utilized. The complexity for solving the main optimization problem is polynomial time, which (regarding the fastest interior-point algorithms) is proportional to K3L, where K is the number of scalar LMI variables and L is the number of scalar LMI rows (see [31]). The complexities of Algorithm 1 in [16], Algorithm 1 in [19], and Algorithm 8 are summarized in Table 1. From Table 1, it can be obtained that the computational burden of the main optimization problem is different mainly in how many of scalar LMI rows are utilized to describe the bounds of estimation error (or true state). In the main optimization problem of [16], L depends on the dimension of true state nx; L in the main optimization problem of [19] can be confined by properly choosing p^. In Algorithm 8, the restriction on the order of zonotope by Nz can limit VX(k).Table 1
The complexity analysis of the compared algorithms.
Different algorithms
Main optimization problem
Refreshment of estimation error set and/or determination of feasibility problem
Algorithm  1 in [16]
K
=
K
M, L=LM+2nx
Each iteration of (35) in [16]: K=nx(nx+1)+1, L=nx2+2nxnpw2nx+2nx(1+nx)+3nx+1; (35) in [16] is infeasible: K=1, L=2nxnpw2nx
Algorithm  1 in [19]
For polyhedral set:K=KM, L=LM+p^;for ellipsoidal set: K=KM+1, L=LM+3
Quickhull Algorithm:O(p^npwlog⁡rp), nx≤3; O(p^npwfr/rp), nx≥4; outer approximation of ellipsoidal sets: K=p^, L=p^2npw
Algorithm8
K
=
K
M, L=LM+VX(k)
Quickhull Algorithm:O(2Nzlog⁡rp), nx≤3; O(2Nzfr/rp), nx≥4; problem (35): K=1, L=2nx+1
Note:KM=(1/2)L(L+1)nx2+Lnx(ny+nu)+nuny+2nx(nx+1)+(1/2)nu(nu+1)+(1/2)ny(ny+1)+1; LM=((L+n-1)!/n!(L-1)!)(4nx+nw+ny+nu)+L((L+n-1)!/n!(L-1)!)(2nx+2nw+ny)+L(2nx+nw+nu)+4nx+nu+ny; npw is the number of vertices related to the outer approximation of εPw by V-presentation; p^ is the maximal number of vertices related to the estimation error set (or true state) in the main optimization problem; Nz is the maximal order of zonotopes related to the estimation error set (or true state) in the main optimization problem.For the auxiliary optimization problem in [16], the computational burden mainly comes from an iterative method. Algorithm 1 in [19] and Algorithm 8 take the advantage of Quickhull Algorithm in [26] to remove the redundant vertices and obtain the estimation error set by V-presentation. The worst-case complexity of Quickhull in [26] is O(Vlog⁡rp) for nx≤3 and O(Vfr/rp) for nx≥4, where V is the number of input points in Rnx, rp is the number of processed points, and fr is the maximum number of facets for rp processed vertices, respectively. Furthermore, the outer approximation of ellipsoidal sets in [19] and problem (35) are solved by an LMI tool, and the complexity analysis of the related parameters K and L is listed in Table 1.Theorem 9.
For system (6), Algorithm 8 is applied. Suppose x(0)∈X(0)=xc(0)⊕H(0)Bm(0), where H(0)∈Rnx×m(0). If the main optimization problem (26) is feasible at time k=0, then x~(k) will converge to a neighborhood of x~=0, and the input and output constraints (8) are satisfied for all time k≥0.Proof.
The proof is similar to the “proof of Theorem 1” in [32]. Suppose, at time k, Algorithm 8 is performed with a refreshed control law; at time k+1,k+2,…,k+l-1,  l>0, it obtains pro=0, and at time k+l it obtains pro=1. According to Algorithm 8, the optimal solutions to (26) at time k will be feasible to (26) at time k+l. Hence, γ(k+l)=γ*(k) is a feasible bound of performance cost at time k+l. At time k+l, solve the main optimization problem (26), γ*(k+l)≤γ(k+l) will be obtained, and then it will satisfy γ*(k+l)≤γ*(k). By induction, γ*(ko) will not increase with time ko, where ko is the time when the main optimization problem (26) is solved. With the evolution of time, γ*(ko) will tend to a constant value denoted by γc, and the corresponding positively invariant ellipsoidal set related to γc for system (11) is represented by ellipsoidal set εQc-1.
Consider the following uncorrupt system (i.e., the system is not corrupted by disturbance):(36)x~ui+1∣ko=Φi,kox~ui∣ko,∀i≥0,x~u0∣ko=x~ko,uui∣ko=Cci∣koxc,ui∣ko+Dci∣koyui∣ko,yui∣ko=Cko+ixui∣ko.The quadratic boundedness condition (23) guarantees the asymptotic stability of the uncorrupt system (36). Hence,(37)x~ui∣koQ-12-x~ui+1∣koQ-12≥1γkoyui∣koQ2+uui∣koR2,kkkkkkkkkkkkkkkikkkkkkkkkkkkkkkk∀i≥0.By summing (37) from i=0 to i=∞, it achieves(38)J∞ko=∑i=0∞yui∣koQ2+uui∣koR2≤γkox~koQ-1≤γko.Hence, γ(ko) is also an upper bound of performance cost for the uncorrupt system (36). The condition J∞(ko)≤γ*(ko) results in limi→∞yu(i∣ko)=0 and  limi→∞uu(i∣ko)=0. When γ*(ko) tends to a constant value, γc means that the augmented state will converge to a neighborhood of x~=0 represented by εQc-1 and stay in this neighborhood thereafter. Satisfaction of the input and output constraints (8) is due to (25).


## 4.1. The Main Optimization Problem

In the dynamic OFRMPC approach, at each timek, based on the dynamic OFRMPC approaches in [16, 19], the following main optimization problem is considered:(15)min⁡γ,Q,A-clj,B-cl,C-cj,D-cmax⁡A∣B∣C∣D∣Ek+i∈Ω,wk+i∈εPwγ,(16)s.t.x~k∈εQ-1,(17)x~i∣kQ-12≥1⟹x~i∣kQ-12-x~i+1∣kQ-12≥1γyi∣kQ2+ui∣kR2,∀i≥0,(18)-u-≤ui∣k≤u-,∀i≥0,(19)-ψ-≤Ψyi+1∣k≤ψ-,∀i≥0,where γ is the performance cost, (16) is the augmented state constraint, (17) is the stability/optimality condition, (18) is the input constraint, (19) is the output constraint, the matrices Q, R are the positive-definite weighting matrices, and Q satisfies the following conditions:(20)Q-1=M=M1M2TM2M3,Q=Q1Q2TQ2Q3.Choosing M2=-M1, then M3=M1(M1-Q1-1)M1, Q2=Q3=Q1-M1-1.In the following, some previous results mainly in [16, 19] are reviewed to deal with constraint conditions (16)–(19). By applying the definition of augmented state and (14), constraint (16) holds if there exist positive-definite matrices {M1,Q1,M-3,M4} such that(21)χrkxckTM1-M1-M1M-3χrkxck≤1,kkkkkkkkkkkkikr∈1,2,…,VXk,(22)M-3-M1-M1M4≥0,M1-M4IIQ1≥0.For all x~(i∣k) and all admissible w(k+i), condition (17) holds if there exist a scalar α and matrices {A^clj,B^cl,C^cj,D^c} such that(23)ΥQBi,k=∑l=1L∑j=1Lλlk+iλjk+iΥljQB≥0,∀i≥0,ΥljQB=1-αM1⋆⋆⋆⋆⋆⋆1-αI1-αQ1⋆⋆⋆⋆⋆00αPw⋆⋆⋆⋆Al+BlD^cCjAlQ1+BlC^cjBlD^cEj+DlQ1⋆⋆⋆M1Al+B^clCjA^cljB^clEj+M1DlIM1⋆⋆Q1/2CjQ1/2CjQ1Q1/2Ej00γI⋆R1/2D^cCjR1/2C^cjR1/2D^cEj000γI,with (10) being parameterized as(24)D-c=D^c,C-cj=D-cCjQ1-C^cjM1-1-Q1-1,B-cl=M1-1M1BlD-c-B^cl,A-clj=M1-1A^clj-AlQ1-BlD-cCjQ1+B-clCjQ1×M1-1-Q1-1+BlC-cj.Suppose, at timek, constrains (21), (22), and (23) are satisfied; there exist positive-definite matrices {U,Ξ~} such that(25)Υju=M1⋆⋆⋆IQ1⋆⋆00Pw⋆2D^cCj2C^cj2D^cEjU≥0,kikUss≤u-s2,j=1,…,L,s=1,…,nu,Υhi,k=∑l=1L∑j=1Lλlk+iλjk+iΥljh≥0,kkkkkkikh=1,…,L,∀i≥0,Υljh=M1⋆⋆⋆IQ1⋆⋆00Pw⋆ΨChAl+BlD^cCjΨChAlQ1+BlC^cjΨChBlD^cEj+DlΞ~,kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkΞ~ss≤ψ~s2,s=1,…,q,ψ~s=12ψ-s-ζs,ζs=max⁡hΨEhPw-1EhTΨTss,where Uss(Ξ~ss) is the sth diagonal element of U(Ξ~), and then by applying (10) and (24), the constraints (18) and (19) are satisfied.From the above conditions, it is shown that the main optimization problem (15)–(19) can be approximated by(26)min⁡γ,α,Q1,M1,M-3,M4,A^clj,B^cl,C^cj,D^c,U,Ξ~γ,is.t.21,22,23,25.At each time k, if the main optimization problem (26) is solved, then (11) is quadratically bounded [29, 30] within a common Lyapunov matrix Q-1, and the controller parameters are calculated by (10) and (24). The main optimization problem (26) can be solved by an LMI (linear matrix inequality) tool, which is the same as in [16, 19].Remark 5.
In [16], the true state constraint in the main optimization problem is represented by the V-representation of general polyhedral sets. In [19], to guarantee the recursive feasibility of the optimization problem based on polyhedral bounds, the true state constraint in the main optimization problem utilizes either the ellipsoidal set or the polyhedral set, with the latter being used if and only if it is contained in the former. Compared with [16, 19], the present paper utilizes the V-representation of zonotopes to represent the bounds of true state.


## 4.2. The Refreshment of Estimation Error Sets

Define the estimation errore(k)≜x(k)-xc(k). At each time k, for the quasi-LPV system (6), the model parameters [A∣B∣C∣D∣E](k)are exactly known at the current time k. If the controller parameters {Ac(k),Bc(k),Cc(k),Dc(k)} are obtained, then according to (11) and the definition of estimation error, it is shown that(27)ek+1=∑l=1L∑j=1Lλlkλjk×Al+BlD-cCj-B-clCjxk+BlC-cj-A-cljxck+BlD-cEj+Dl-B-clEjwk=Ak+BkDckCk-BckCkxk+BkCck-Ackxck+BkDckEk+Dk-BckEkwk.By applying (7) and (13) and Properties 1 and 2, it can be obtained that(28)ek+1∈Ek+1≜Ak+BkDckCk-BckCk+BkCck-Ackxck⊕Ak+BkDckCk-BckCkHk,BkDckEk+Dk-BckEkHw×Bmk+d.From (28), the estimation error set at time k+1 can be represented by a [m(k)+d]-zonotope E(k+1)=p(k+1)⊕H-(k+1)Bm(k)+d, where(29)pk+1=Ak+BkDckCk-BckCk+BkCck-Ackxck,H-k+1=Ak+BkDckCk-BckCkHk,BkDckEk+Dk-BckEkHw.It can easily be obtained from (28) that, with the evolution of time, the order of zonotope E(k+1) will increase. Accordingly, the generators and the vertices of zonotope E(k+1) will increase, which may lead to more LMIs in the main optimization problem (26) to enforce the augmented state constraints and then increase the online computational burden of the main optimization problem (26). To decrease the online computational burden of the main optimization problem (26), we impose the restriction on the order of E(k+1) by m(k)+d≤Nz. When m(k)+d>Nz, Property 3 is applied to construct a low order Nz-zonotope E^(k+1)=p(k+1)⊕H^(k+1)BNz for outer approximation of the high order zonotope E(k+1); that is, E(k+1)⊆E^(k+1), where(30)H^k+1=H^Tk+1Ok+1,kkikkH^k+1∈Rnx×Nz,and H^T(k+1)∈Rnx×[Nz-nx] and O(k+1)∈Rnx×nx are calculated according to (2) and (3), respectively. From the above derivations, it can be obtained that the estimation error set at time k+1 is one of the following forms:(31)Pk+1=Ek+1,Hk+1=H-k+1,mk+1=mk+d,whenmk+d≤Nz;Pk+1=E^k+1,Hk+1=H^k+1,mk+1=Nz,whenmk+d>Nz.When the estimation error setP(k+1) is obtained, according to Definition 1, at most 2m(k+1) vertices can be calculated. By removing the redundant vertices, the set P(k+1) can be described by V-representation through VP(k+1) vertices. The vertices of the true state set X(k+1) are the same as P(k+1); that is, VX(k+1)=VP(k+1). Hence the V-representation of the true state at time k+1 is(32)xk+1∈Coχ1k+1,χ2k+1,…,χVXk+1k+1,χrk+1=xck+1+vrk+1,kkir∈1,2,…,VXk+1.Remark 6.
In Section2, we have mentioned the relationship between Vz (the number of vertices related to zonotopes), m (the order of zonotopes), and the dimension of zonotopes. For the refreshment of estimation error sets in Section 4.2, when the dimension of true state is nx=2, VX(k+1)=VP(k+1)=2m(k+1); when nx=3, VX(k+1)=VP(k+1)=m(k+1)[m(k+1)-1]+2, respectively [27]. For high-dimensional systems, when nx is fixed, the upper bound of VX(k+1) (or VP(k+1)) is related to m(k+1) and can be calculated in polynomial time [28]. In [16], VX(k+1) (or VP(k+1)) only depends on the dimension of true state; that is, VX(k+1)=VP(k+1)=2nx. In [19], the upper bound of VX(k+1) (or VP(k+1)) is bounded by properly selecting an integer p^.
In [16, 19], the outer approximation of ellipsoidal bounded disturbance is represented by the V-representation of the polytope. The present paper utilizes zonotopes to outer approximation of εPw. The outer approximation of ellipsoids introduces conservatism. However, by properly increasing the vertices of outer approximation polyhedral sets in [16, 19] and choosing matrix Hw with proper generators, the conservatism can be reduced.


## 4.3. The Feasibility of the Optimization Problem at the Next Sampling Time

The original optimization problem (15)–(19) is recursive feasible; that is, it is feasible for all k>0 if and only if it is feasible at time k=0. However, by using the bounds of true state to handle the uncertainty of true state, the recursive feasibility may be lost for the main optimization problem (26). If the main optimization problem (26) is solved at time k, then the optimal solutions {α,γ,A^clj,B^cl,C^cj,D^c,Q1,M1,U,Ξ~}*(k) are obtained, and x~(k)∈ε[Q*(k)]-1. At time k+1, by choosing(33)M4k+1=M1*k-Q1*k-1,M1k+1=M1*k,Q1k+1=Q1*k,M-3k+1=M3k+1=M1*kM1*k-Q1*k-1-1M1*kT,then (22) is satisfied as follows:(34)M-3k+1-M1k+1-M1k+1M4k+1≥0,M1k+1-M4k+1IIQ1k+1≥0.Based on (33) and (34) and the definitions of Q and M, then, by applying the Schur complement, it can be obtained that if the problem (35) is feasible:(35)findη,s.t.η⋆⋆χrk+1Q1*k⋆xck+1Q2*kQ3*k≥0,mmi0<η≤1,r∈1,2…,VXk+1,then (21) is satisfied and x~(k+1)∈ε[Q*(k)]-1. In the main optimization problem (26), only (21) and (22) are related to the process state information. By further choosing {α,γ,A^clj,B^cl,C^cj,D^c,U,Ξ~}(k+1)={α,γ,A^clj,B^cl,C^cj,D^c,U,Ξ~}*(k), the feasible solutions for the main optimization problem (26) at timek+1can be constructed. The feasible solutions for (21) at timek+1can be directly constructed through (33). Hence, the feasibility of (35) determines not only whether x~(k+1) is contained in the ellipsoidal set ε[Q*(k)]-1 or not, but also the feasibility of the main optimization problem (26) at time k+1.Remark 7.
In [16], the refreshment of estimation error sets and the feasibility of the main optimization problem at the next sampling time are performed simultaneously by an iterative auxiliary optimization problem. Furthermore, it involves the H-representation of estimation error set in the auxiliary optimization problem. In [19], to guarantee the recursive feasibility of the main optimization problem based on polyhedral sets, the recursions of both polyhedral sets and ellipsoidal sets are obtained. At the next sampling time, the polyhedral set is used in the main optimization problem if and only if it is contained in the ellipsoidal set. The present paper refreshes the estimation error set via zonotopic computations and is represented by V-representation. Then problem (35) is applied to determine whether the main optimization problem (26) is solved or not at the next sampling time.


## 4.4. Overall Solutions

Algorithm 8.
Choose a feasiblexc(0) and the initial constrain of x(0) represented by a zonotope xc(0)⊕H(0)Bm(0), where H(0)∈Rnx×m(0). Let pro=1. Select the maximal order of zonotope related to the estimation error set as Nz. At each time k≥0, perform the following steps.(1)
Ifpro=1, then solve the main optimization problem (26) and obtain the controller parameters by (10) and (24).(2)
Ifpro=0, then inherit the previous controller parameters.(3)
Implementu(k)=Cc(k)xc(k)+Dc(k)y(k), and calculate xc(k+1)=Ac(k)xc(k)+Bc(k)y(k).(4)
CalculateE(k+1) by (28). If m(k)+d≤Nz, then let P(k+1)=E(k+1); else then calculate E^(k+1) by (30) and let P(k+1)=E^(k+1).(5)
Check the feasibility of problem (35). If problem (35) is feasible, then letpro=1; else then let pro=0.The complexity analysis of Algorithm8 and the related algorithms for worst cases are listed in Table 1. For the main optimization problem in Table 1, an LMI tool is utilized. The complexity for solving the main optimization problem is polynomial time, which (regarding the fastest interior-point algorithms) is proportional to K3L, where K is the number of scalar LMI variables and L is the number of scalar LMI rows (see [31]). The complexities of Algorithm 1 in [16], Algorithm 1 in [19], and Algorithm 8 are summarized in Table 1. From Table 1, it can be obtained that the computational burden of the main optimization problem is different mainly in how many of scalar LMI rows are utilized to describe the bounds of estimation error (or true state). In the main optimization problem of [16], L depends on the dimension of true state nx; L in the main optimization problem of [19] can be confined by properly choosing p^. In Algorithm 8, the restriction on the order of zonotope by Nz can limit VX(k).Table 1
The complexity analysis of the compared algorithms.
Different algorithms
Main optimization problem
Refreshment of estimation error set and/or determination of feasibility problem
Algorithm  1 in [16]
K
=
K
M, L=LM+2nx
Each iteration of (35) in [16]: K=nx(nx+1)+1, L=nx2+2nxnpw2nx+2nx(1+nx)+3nx+1; (35) in [16] is infeasible: K=1, L=2nxnpw2nx
Algorithm  1 in [19]
For polyhedral set:K=KM, L=LM+p^;for ellipsoidal set: K=KM+1, L=LM+3
Quickhull Algorithm:O(p^npwlog⁡rp), nx≤3; O(p^npwfr/rp), nx≥4; outer approximation of ellipsoidal sets: K=p^, L=p^2npw
Algorithm8
K
=
K
M, L=LM+VX(k)
Quickhull Algorithm:O(2Nzlog⁡rp), nx≤3; O(2Nzfr/rp), nx≥4; problem (35): K=1, L=2nx+1
Note:KM=(1/2)L(L+1)nx2+Lnx(ny+nu)+nuny+2nx(nx+1)+(1/2)nu(nu+1)+(1/2)ny(ny+1)+1; LM=((L+n-1)!/n!(L-1)!)(4nx+nw+ny+nu)+L((L+n-1)!/n!(L-1)!)(2nx+2nw+ny)+L(2nx+nw+nu)+4nx+nu+ny; npw is the number of vertices related to the outer approximation of εPw by V-presentation; p^ is the maximal number of vertices related to the estimation error set (or true state) in the main optimization problem; Nz is the maximal order of zonotopes related to the estimation error set (or true state) in the main optimization problem.For the auxiliary optimization problem in [16], the computational burden mainly comes from an iterative method. Algorithm 1 in [19] and Algorithm 8 take the advantage of Quickhull Algorithm in [26] to remove the redundant vertices and obtain the estimation error set by V-presentation. The worst-case complexity of Quickhull in [26] is O(Vlog⁡rp) for nx≤3 and O(Vfr/rp) for nx≥4, where V is the number of input points in Rnx, rp is the number of processed points, and fr is the maximum number of facets for rp processed vertices, respectively. Furthermore, the outer approximation of ellipsoidal sets in [19] and problem (35) are solved by an LMI tool, and the complexity analysis of the related parameters K and L is listed in Table 1.Theorem 9.
For system (6), Algorithm 8 is applied. Suppose x(0)∈X(0)=xc(0)⊕H(0)Bm(0), where H(0)∈Rnx×m(0). If the main optimization problem (26) is feasible at time k=0, then x~(k) will converge to a neighborhood of x~=0, and the input and output constraints (8) are satisfied for all time k≥0.Proof.
The proof is similar to the “proof of Theorem 1” in [32]. Suppose, at time k, Algorithm 8 is performed with a refreshed control law; at time k+1,k+2,…,k+l-1,  l>0, it obtains pro=0, and at time k+l it obtains pro=1. According to Algorithm 8, the optimal solutions to (26) at time k will be feasible to (26) at time k+l. Hence, γ(k+l)=γ*(k) is a feasible bound of performance cost at time k+l. At time k+l, solve the main optimization problem (26), γ*(k+l)≤γ(k+l) will be obtained, and then it will satisfy γ*(k+l)≤γ*(k). By induction, γ*(ko) will not increase with time ko, where ko is the time when the main optimization problem (26) is solved. With the evolution of time, γ*(ko) will tend to a constant value denoted by γc, and the corresponding positively invariant ellipsoidal set related to γc for system (11) is represented by ellipsoidal set εQc-1.
Consider the following uncorrupt system (i.e., the system is not corrupted by disturbance):(36)x~ui+1∣ko=Φi,kox~ui∣ko,∀i≥0,x~u0∣ko=x~ko,uui∣ko=Cci∣koxc,ui∣ko+Dci∣koyui∣ko,yui∣ko=Cko+ixui∣ko.The quadratic boundedness condition (23) guarantees the asymptotic stability of the uncorrupt system (36). Hence,(37)x~ui∣koQ-12-x~ui+1∣koQ-12≥1γkoyui∣koQ2+uui∣koR2,kkkkkkkkkkkkkkkikkkkkkkkkkkkkkkk∀i≥0.By summing (37) from i=0 to i=∞, it achieves(38)J∞ko=∑i=0∞yui∣koQ2+uui∣koR2≤γkox~koQ-1≤γko.Hence, γ(ko) is also an upper bound of performance cost for the uncorrupt system (36). The condition J∞(ko)≤γ*(ko) results in limi→∞yu(i∣ko)=0 and  limi→∞uu(i∣ko)=0. When γ*(ko) tends to a constant value, γc means that the augmented state will converge to a neighborhood of x~=0 represented by εQc-1 and stay in this neighborhood thereafter. Satisfaction of the input and output constraints (8) is due to (25).


## 5. Numerical Example

Consider the nonlinear continuous stirred tank reactor (CSTR) model, which has been studied in [12, 16, 19]. Similar to [16, 19], the nonlinear CSTR model is represented as the quasi-LPV system around the operation point, which considers the physical constraints and bounded disturbance. The details corresponding to system (6) are L=4,(39)A1=0.8227-0.001686.12330.9367,A2=0.9654-0.00182-0.67590.9433,A3=0.8895-0.002942.94470.9968,A4=0.8930-0.000622.77380.8864,B1=-0.0000920.1014,B2=-0.0000970.1016,B3=-0.0001570.1045,B4=-0.0000340.0986,Dl=0.00220.0564,Cl=01,El=0.2,kkkkkkkkkkkkkkkkkkkkkkkkkkkkl=1,…,4,u-=10,Ψ=1,ψ-=10,λ1=12φ1y-φ1-ψ-φ1ψ--φ1-ψ-,λ2=12φ1ψ--φ1yφ1ψ--φ1-ψ-,kkkkkkkkkkkkkkkkkkφ1y=7.2×1010e-8750/y+350,λ3=12φ2y-φ2-ψ-φ2ψ--φ2-ψ-,λ4=12φ2ψ--φ2yφ2ψ--φ2-ψ-,kkkkkkkφ2y=3.6×1010e-8750/y+350-e-8750/350y.For the above quasi-LPV system, in Algorithm8, choose Pw=4, Hw=0.5, Q=25, and R=1, take xc(0)=[0.15,3]T, and the initial state belongs to a zonotope xc(0)⊕H(0)Bm(0), where H(0)=0.12001.2, m(0)=2. In order to compare with Algorithm 1 in [16], Nz=2 (VX(k)=VP(k)=4) is selected in Algorithm 8, and choose the parameters for Algorithm 1 in [16] as e-=[0.12,1.2]T, {τ-,δ,μ}={10,0.003,0.99}. For comparing with Algorithm 1 in [19], Nz=4 (VX(k)=VP(k)=8) is selected in Algorithm 8, and choose the parameters for Algorithm 1 in [19] as e-=[0.12,1.2]T, T=0.16970.16970-1.69700.16971.6971.697T, Me(0)=diag{0.5/(0.06×2)2,0.5/(0.6×2)2},  p^=8.For all the compared algorithms above, the disturbance sequencew is randomly generated from the interval [-0.5,0.5], which can be represented by a zonotope D=0⊕HwB1. We have fixed α=0.001, which only has negligible influence on the control performance. By applying Algorithm 8 (Nz=2), Algorithm 8 (Nz=4), Algorithm 1 in [16], and Algorithm 1 in [19], the simulation results are shown in Figures 2, 3, 4, 5, 6, 7, 8, 9, 10, and 11. Figures 2–4 show the responses of estimated states and time-varying estimation error sets related to the compared algorithms. Figure 6 compares the trajectories of augmented closed-loop system corresponding to the different algorithms. The control input signals shown in Figure 7 illustrate that the input constraints are satisfied. Figure 8 shows the performance cost of the compared algorithms. The same disturbance sequence that is shown in Figure 9 has been applied for all the simulations. Figures 10 and 11 show the projections of a super dimensional ellipsoid εQc-1, which is related to performance cost γc, on the subspace of estimated state and estimation error set, respectively.Figure 2
The responses of estimated states and time-varying estimation error sets, Algorithm 1 in [16].Figure 3
The responses of estimated states and time-varying estimation error sets, Algorithm 1 in [19].Figure 4
The responses of estimated states and time-varying estimation error sets, Algorithm8, Nz=2.Figure 5
The responses of estimated states and time-varying estimation error sets, Algorithm8, Nz=4.Figure 6
The comparisons of trajectories of augmented closed-loop systems.Figure 7
The comparisons of control input signals.Figure 8
The performance costs of the compared algorithms.Figure 9
The disturbance utilized in the simulations.Figure 10
The projection ofεQc-1 on the subspace of estimated state.Figure 11
The projection ofεQc-1 on the subspace of estimation error.LetJ=∑k=030{y(k)Q2+u(k)R2}. For Algorithm 8 (Nz=2) (Algorithm 8 (Nz=4), Algorithm 1 in [16], and Algorithm 1 in [19]), J=3808.8 (J=3374.5, J=4253.6, J=3658.8) is obtained, and the time spent on the simulations is 36.4 s (42.6s,50.8s,49.2s), respectively. Compare Algorithm 8 (Nz=2) with Algorithm 1 in [16]; the control performance is improved with lower computational burden; even the main optimization problem of Algorithm 1 in [16] is not solved at some sampling times. Compare Algorithm 8 (Nz=2) with Algorithm 8 (Nz=4); by properly limiting the order of zonotopes, the control performance can be improved, with a slight increasement of computational burden. Furthermore, the augmented states of Algorithm 8 (Nz=4) and Algorithm 1 in [19] converge to nearly the same neighborhoods of x~=0, while the control performance of the former is better than the latter. The Matlab 7.14 (AMD Phenom II Processor 2.70 GHz, 4 G Memory) is utilized for the simulations.


## 6. Conclusions

For quasi-LPV systems with bounded disturbance, a set-membership estimation method via zonotopic computation is applied to refresh the estimation error set. By properly refreshing the estimation error set, the control performance can be improved. Furthermore, the feasibility of the main optimization problem at the next sampling time can be determined by checking a feasible problem at the current time. The augmented state can converge to a neighborhood ofx~=0 and stay in the neighborhood of x~=0 thereafter.



---

*Source: 875850-2015-03-03.xml*