A Robust Constrained Reference Governor Approach using Linear Matrix Inequalities

A Robust Constrained Reference Governor Approach using Linear Matrix Inequalities J.L. Guzm´an a∗,T.´ Alamo b, M. Berenguel a,S.Dormidoc, E. F. Camacho b aDep. Lenguajes y Computaci´on, Universidad de Almer´ıa, Ctra. Sacramento s/n 04120 Almer´ıa, Spain. E-mail:jo[email protected] bDep. Ingenier´ıa de Sistemas y Autom´atica, Universidad de Sevilla, C. de los Descubrimientos s/n 41092 Sevilla, Spain. cDep. Inform´atica y Autom´atica, ETSI Inform´atica, UNED, 28040 Madrid, Spain. Abstract The purpose of this paper is to examine and provide a solution to the output reference tracking problem for uncertain systems subject to input saturation. As well-known, input saturation and modelling errors are very common problems at industry, where control schemes are implemented without accounting for such problems. In many cases, it is sometimes difficult to modify the existing implemented control schemes being necessary to provide them with external supervisory control approaches in order to tackle problems with constraints and modelling errors. In this way, a cascade structure is proposed, combining an inner loop containing any controller with an outer loop where a Generalized Predictive Controller (GPC) provides adequate references for the inner loop considering input saturations and uncertainties. Therefore, the contribution of this paper consists in providing a state space representation for the inner loop and using Linear Matrix Inequalities (LMI) to obtain a predictive state-vector feedback in such a way that the input reference for the inner loop is calculated to satisfy robust tracking specifications considering input saturations. Hence, the final proposed solution consists in solving a regulation problem to a fixed reference value subjected to a set of constraints described by several LMI and Bilinear Matrix Inequalities (BMI). An illustrative numerical example is presented. Key words: predictive control, robust control, tracking control, constrained control, linear matrix inequalities 1. Introduction Mathematical models are required at design level in any control system development. Models cannot represent every aspect of reality, so assumptions must be made in order to use them for control purposes. On the other hand, most physical processes are constrained by several reasons, such as physical limits (e.g. valve position), security levels (e.g. pressure levels), or performance criteria (e.g. working near the optimal operating point) [1]. In practice, many control techniques implemented at industry work without taking into account these modelling errors and system constraints. Fixed-structure models and known parameters are used supposing that the model exactly represents the real process, and the imperfections will be removed by means of feedback. Furthermore, detuned control is usually used in order to cope with system saturations. As well-know, these problems have widely been addressed by the scientific control community. Numerous robust control techniques are available in order to face system uncertainties such as H∞[2], Quantitative Feedback Theory (QFT) [3], or μ-synthesis [4]. In the same way, constrained system problems have been solved from different points of view including anti-windup schemes [5],[6], constrained model predictive control [7], [8], and LMI-based synthesis [9]. Furthermore, combinations of these techniques can be found in order to solve both robustness and input saturation problems [10], [11], [12], [13]. All previous approaches can be used in order to control systems presenting the problems described above. However, many industrial processes are currently controlled by some traditional control schemes such as PID control or a generic two degrees of freedom controller, being difficult or even impossible to modify this primal controller. Therefore, in these cases, an external supervisory control is required in order to face the problems of the implemented inner loop. In this sense, the reference governor approach has been presented in several works as solution to this problem, mainly focused on providing appropriated references to the inner loop in order to ensure nominal stability in presence of constraints [14], [15], [16], [17], [18]. This approach has also been studied for solving robustness problems due to the presence of uncertainties in the inner loop [19], [20], Preprint submitted to Journal of Process Control 17 April 2007 Manuscript [21]. However, the existing reference governor approaches present some drawbacks, such as the requirement of knowing the future reference in advance, the use of very conservative uncertainty representations, the heavy computational burden or the difficulty to improve tracking performance. Therefore, this work presents a robust constrained reference governor approach based on LMI with the aim of solving the previous problems in a systematic way [22]. In this approach a GPC controller [23], [7] provides adequate references for the inner loop considering input saturations and uncertainties. Thus, the proposed solution in order to prove robust stability when constraints are active, is to translate the inner loop problem into a state space representation, and then using LMI to obtain a predictive state-vector feedback in such a way that the input reference to this inner loop is calculated in order to satisfy robust tracking specifications considering input saturation. This work describes the different steps required to obtain such solution and how input saturation in the inner loop can be handled using LMI-based methods. As it will be shown, non-symmetric limits are also allowed. The final solution consists in solving a set of constraints defined by several LMI and BMI, where a Branch and Bound algorithm has been developed in order to handle the bilinear terms. It is important to notice that the algorithm is implemented for tracking problems where the aim is to regulate to a fixed reference value and not to the origin, in the presence of input constraints in the inner loop of the system what is not usual in the referenced works. The paper is organized as follows. Section 2 describes the state space representation of the inner loop including the saturation term. The next section presents the LMI-based solution for this approach ensuring constrained robust stability. The problem is formulated by means of LMI in order to ensure constrained robust stability and fulfill specific performance criteria, describing also the proposed Branch and Bound algorithm to handle the bilinear terms. Finally, a numerical example, which can be found in many industrial plants, is presented in Section 4. 2. State space representation of the inner loop. The first step necessary before applying a LMI-based solution is to obtain the state space representation for the inner loop in the proposed approach (see Figure 1). The inner loop has been considered as a typical control scheme with two degrees of freedom for general purposes. In order to provide optimal references to the inner loop, the GPC algorithm has been selected. GPC is based on CARIMA model, where the following plant representation is considered being the time delay included into the B(z−1) polynomial A(z−1)y(t)=B(z−1)u(t−1) (1) with A(z−1)= n  i=1 (1 + piz−1),B(z−1)=Kz m  i=1 (1 + ciz−1) (2) Fig. 1. Control system scheme for LMI-based approach The plant parameters Kz,c1,...,c m,p1,...,p nare considered uncertain and supposed to lie within known intervals due to the presence of uncertainties in the plant time domain parameters. That is, Kz∈[Kz,min,K z,max] ci∈[ci,min,c i,max],i=1,...,m pi∈[pi,min,p i,max],i=1,...,n The polynomials A(z−1)andB(z−1) can be rewritten as A(z−1)=1+a1z−1+a2z−1+...+anz−n(3) B(z−1)=b0+b1z−1+b2z−1+...+bmz−m(4) where the coefficients a1,...,a n,b0,...,b m, depend multilineally on the parametric vector φ=[Kz,c 1,...,c m,p 1,...,p n]T . That is, each coefficient depends affinely on each element of vector φ. In order to express explicitly the dependence on the vector φ, the coefficients can be expressed as b0(φ),...,b m(φ), a1(φ),...,a n(φ). Considering the plant dynamics (1) and the polynomials A(z−1)andB(z−1), it results that y(t)=− n  i=1 ai(φ)y(t−i)+ m  i=0 bi(φ)u(t−i−1) (5) or y(t+1)=− n  i=1 ai(φ)y(t−i+1)+ m  i=0 bi(φ)u(t−i) (6) The plant dynamics can be represented by a state-space representation, where the proposed state depends on the current output, and the past outputs and inputs in the following way xp(t)=[y(t)... y(t−n+1)u(t−1) ... u(t−m)]T(7) This state selection has the advantage that the state xp(t)is always accessible, that is, the value of xp(t) is known since it is always possible to access to the output y(t) and input u(t) signals. So, the state space representation is given by x+ p=Apxp+Bpu(8) 2 y=Cpxp(9) where Ap= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ −a1−a2... −an−1−anb1b2... bm−1bm 10... 0000... 00 01... 00 00... 00 . . .. . ..... . .. . .. . .. . ..... . .. . . 00... 1000... 00 00... 0000... 00 00... 00 10... 00 00... 00 01... 00 . . .. . ..... . .. . .. . .. . ..... . .. . . 00... 0000... 10 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,Bp= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ b0 0 . . . 0 0 1 0 0 . . . 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Cp=100... 0 xpdenotes the state vector, uthe system input, and x+ pthe next state xp(t+ 1). Considering that the coefficients a1,...,a n,b0,...,b m depend on the parametric vector φ, the system dynamics can be rewritten as x+ p=Ap(φ)xp+Bp(φ)u(10) y=Cpxp(11) The plant disturbances have not been included along the previous description. Bounded disturbances in the plant represented by the CARIMA model are given as A(z−1)y(t)=B(z−1)u(t−1) + T(z−1)(t) Δ(12) where Δ = 1 −z−1. It is assumed that the coefficients of the T(z−1) polynomial depends multilineally on a bounded parametric vector φTand (t) is bounded for all t>0, that is, (t)∞< max,∀t. Following the same steps as previously, the system dynamics can be rewritten as x+ p=Ap(˜ φ)xp+Bp(˜ φ)u+Ep(˜ φ) (13) y=Cpxp In this representation, ˜ φis a parametric vector containing φ,φTand . Furthermore, it can be assumed that ˜ φ can only take values within a convex set (typically an hyperrectangle). Finally, notice that Ap(˜ φ), Bp(˜ φ)andEp(˜ φ) depends multilineally on parametric vector ˜ φ. 2.0.1. Controller and prefilter representation. Assume available state space descriptions for the prefilter F(z−1) and controller C(z−1). Denoting xFas the state vector of the filter F(z−1), rthe filter input and rFthe filter output, it is supposed that matrices AF,BF,CF,and DFdescribe the filter dynamics as follows x+ F=AFxF+BFr(14) rF=CFxF+DFr In the same way, xCdenotes the state vector for the controller C(z−1)anduthe controller output. The matrices AC,BC,CCand DCdescribe the controller dynamics as follows x+ C=ACxC+BC(rF−y) (15) u=CCxC+DC(rF−y) Note that the input to the controller is given by the filter output rFminus the plant output y, and the plant is subjected to uncertainties and disturbances as discussed above. 2.0.2. Inner loop representation. As commented previously, the goal is to design a robust predictive controller considering input saturation in the inner loop. Therefore, the state space representation of the inner loop must be developed including the saturation. The input saturation in the inner loop is given by σp(u)=⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Umin if u<U min uif Umin ≤u≤Umax Umax if u>U max (16) where nonsymmetric saturation can be present. Firstly, the saturation is redefined in order to use a symmetric representation to facilitate the calculations. Therefore, the saturation is obtained as σp(u)=Lsσ(1 Ls (u−uc)) + uc(17) where σ(u)=−1if u < −1 uif−1≤u≤1 1if u > 1 (18) uc=Umax +Umin 2,L s=Umax −Umin 2 Then, the plant representation (13) is modified to consider input saturation in the following way x+ p=Apxp+Bp(Lsσ(1 Ls (u−uc)) + uc)+Ep(19) y=Cpxp where Ap=Ap(˜ φ), Bp=Bp(˜ φ), and Ep=Ep(˜ φ) will be considered from now on for the sake of simplicity. The proposed extended vector xincluding the inner loop dynamics is defined as 3 x=⎡ ⎢ ⎢ ⎢ ⎣ xp xC xF ⎤ ⎥ ⎥ ⎥ ⎦ (20) Then, the full system described by the plant, prefilter, and controller has ras input (prefilter input), and yas output (plant output). Then, x+ Fand x+ Ccan be described as function of xand ras follows x+ F=AFxF+BFr=00AFx+BFr(21) x+ C=ACxC+BC(rF−y) = (22) =ACxC+BC(CFxF+DFr)−BCCpxp= =−BCCPxp+ACxC+BCCFxF+BCDFr= =−BCCpACBCCFx+BCDFr being the control signal uobtained in the following way u=CCxC+DC(rF−y) = (23) =CCxC+DC(CFxF+DFr)−DCCpxp= =−DCCpxp+CCxC+DCCFxF+DCDFr= =−DCCpCCDCCFx+DCDFr= =Cux+Dur In this way, using the new plant representation (19), and the prefilter and controller state-space dynamics (21-22), the closed-loop state space representation for the inner loop is described as x+=Ax +BuσCu Ls x+Dur−uc Ls+E+Brr(24) y=Cyx where A=⎡ ⎢ ⎢ ⎢ ⎣ Ap00 −BCCpACBCCF 00AF ⎤ ⎥ ⎥ ⎥ ⎦ ,B u=⎡ ⎢ ⎢ ⎢ ⎣ LsBp 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ , E=⎡ ⎢ ⎢ ⎢ ⎣ Ep+Bpuc 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ ,B r=⎡ ⎢ ⎢ ⎢ ⎣ 0 BCDF BF ⎤ ⎥ ⎥ ⎥ ⎦ ,C y=100  3. Robust Constrained LMI-based approach. 3.1. Tracking problem. Preliminary ideas. Most of the results obtained for constrained MPC using LMI have been proposed to regulate the system to the origin. In this way, the results obtained in [24] can be used to calculate a control law r=Ksxfor the plant (24) considering the system free of disturbances (E(˜ φ) = 0) and regulating to the origin. Notice that x= 0 is an equilibrium for the system and for all value of ˜ φ. Therefore, Kscan be calculated ensuring robust stability and in such a way that the control law r=Ksxregulates to the origin for all posible initial conditions and any value of ˜ φ[24], [25]. However, one of the objectives considered in the approach presented in this work is to make the output yreach the reference value w. Therefore, the problem formulation must be oriented to this objective. This work presents some preliminary ideas based on the extensions proposed in [24] for set-point tracking. Firstly, it is necessary to notice that due to the dependence on the parametric vector ˜ φ, it is imposible to find static values for xand r(xeand re) such that the system finds a unique equilibrium for all values of the parametric vector ˜ φ. In this work in order to address this problem, the following control law is proposed r=re+Ks(x−xe) (25) where xe,re,andKswill be obtained in such a way that the performance of the closed-loop system is enhanced and the system evolution is ensured to be inside an invariant ellipsoid containing the problem initial conditions. Substituting rin equation (24) by the desired control law r=re+Ks(x−xe), the following expression is obtained x+=Ax +BuσCux+Du(re+Ks(x−xe))−uc Ls+E+ (26) +Br(re+Ks(x−xe)) Then, if the change ¯x=x−xeis considered, and so ¯x+=x+−xe, the system dynamics can be represented as ¯x+=A(¯x+xe)+Buσ(Cu+DuKs)¯x+Cuxe+Dure−uc Ls+E+ +Br(re+Ks¯x)−xe y=Cyx=Cy¯x+Cyxe Finally, if the following changes are performed Ku=Cu+DuKs,d u=Cuxe+Dure−uc it is obtained that ¯x+=ABc¯x+BBcσKu¯x+du Ls+EBc (27) y=Cy¯x+Cyxe where ABc =(A+BrKs),B Bc =Bu,E Bc =E+Axe+Brre−xe Define ℘(Ps,ρ)=¯xTPs¯x≤ρas an ellipsoid where ¯x0∈ ℘(Ps,ρ) with ¯x0=x0−xe. In this way, the tracking problem assuring constrained robust stability for the system (27) will be solved using LMI and fulfilling the following objectives: (i) Firstly, the decision variables xe,re,Ks,Ps,andρare calculated in such a way that the ellipsoid ℘(Ps,ρ)is 4 invariant containing the system initial conditions ¯x0 and using the control law r=re+Ks(x−xe). (ii) After that, new constrains will be included in order to fulfill a certain performance criteria. These objectives will be addressed in next sections, but the way in which the saturation term presented in (27) can be taken into account will be addressed before. 3.1.1. Linear Difference Inclusion of the saturation term Notice that due to the input saturation, a nonlinear term appears in the system dynamics, σ(Ku¯x+du Ls). This nonlinear term can be approximated using the Linear Difference Inclusion (LDI) results obtained in [26] and [27] where it is shown that, if b∈Rsatisfies |b|≤1 then σ(a)∈Co{a, b},∀a∈R being Co the convex hull. In particular, if |Hs¯x+h|≤1 ∀¯x∈℘(Ps,ρ) then σ(Ku¯x+du Ls)∈Co{Ku¯x+du Ls,H s¯x+h},∀(Ku¯x+du Ls)∈R,∀¯x∈℘(Ps,ρ) Therefore, each objective commented above will be translated to analyze if it satisfies the extremes of the convex hull ¯x+=ABc¯x+BBcKu¯x+du Ls+EBc (28) ¯x+=ABc¯x+BBc(Hs¯x+h)+EBc (29) as will be shown in next sections. On the other hand, the inequality |Hs¯x+h|≤1 must be considered. This inequality can be translated to a LMI in order to be included in the final optimization problem. The inequality can be expressed as two inequalities in the following way Hs¯x+h≤1⇒Hs¯x≤1−h, ∀¯x∈℘(Ps,ρ) (30) Hs¯x+h≥−1⇒Hs¯x≥−1−h, ∀¯x∈℘(Ps,ρ) (31) where these inequalities must be satisfied in the ellipsoid ℘(Ps,ρ). In the next section, this ellipsoid will be forced to be invariant containing the system initial conditions. Considering the first inequality (30) and using the S− procedure, it is equivalent to study the existence of λ2≥0 such that (Farkas lemma [25]) Hs¯x+h+λ2(ρ−¯xTPs¯x)≤1,∀¯x(32) This can be expressed as ⎡ ⎣ 1 −¯x⎤ ⎦ T⎡ ⎢ ⎣ λ2ρ+h−1−1 2Hs −1 2HT s−λ2Ps ⎤ ⎥ ⎦⎡ ⎣ 1 −¯x⎤ ⎦≤0,∀¯x(33) or equivalently ⎡ ⎢ ⎣ 1−h−λ2ρ1 2Hs 1 2HT sλ2Ps ⎤ ⎥ ⎦>0 (34) Then, pre- and post-multiplying by diag[I2P−1 s] and making W=P−1 s,V=HsW,itresults ⎡ ⎣ 1−h−λ2ρV VT4λ2W⎤ ⎦>0 (35) The same procedure can be applied for the second inequality (31) obtaining ⎡ ⎣ 1+h−λ2ρV VT4λ2W⎤ ⎦>0 (36) For a fixed value of λ2, notice that the obtained matrix inequalities are LMI in the decision variables ρ,Ps,Hs, and h. The procedure to find correct values for λ2will be described later. As consequence of the previous results, the following property is proposed: Property 1 Suppose that there exits λ2≥0such that LMI (35) and (36) are fulfilled, then: σ(Ku¯x+du Ls )∈Co{Ku¯x+du Ls ,H s¯x+h},∀¯x∈℘(Ps,ρ), where P=W−1and H=W−1V. 3.2. Robust invariant ellipsoid As commented previously, one of the objectives is to calculate the decision variables xe,re,Ks,Ps,andρin such a way that the ellipsoid ℘(Ps,ρ) is invariant including the system initial conditions and using the control law r= re+Ks(x−xe). Therefore, in order to ensure the ellipsoid being invariant the following inequality must be fulfilled (¯x+)TPs(¯x+)≤ρ, ∀¯x∈℘(Ps,ρ) (37) This problem can be reformulated using S−procedure as follows: (¯x+)TPs(¯x+)+λ1(ρ−¯xTPs¯x)≤ρ, ∀¯x, λ1≥0 (38) or equivalently (¯x+)TPs(¯x+)−λ1¯xTPs¯x+ρ(λ1−1) ≤0,∀¯x, λ1≥0 (39) Lets consider the following property: Property 2 Suppose that Ps>0, then ϑTPsϑ≥vTPsv+2vTPs(ϑ−v)=−vTPsv+2vTPsϑ and ϑTPsϑ= max v{−vTPsv+2vTPsϑ}(40) Therefore, using the previous property and the closed-loop system dynamics (27), the inequality (39) results −vTPsv+2vTPs(ABc¯x+BBcσ(Ku¯x+du Ls)+EBc)−(41) −λ1¯xTPs¯x+ρ(λ1−1) ≤0 where this inequality must be satisfied ∀¯xand ∀v. 5 Notice that in order to address this problem and demonstrate that the system evolution belongs to an invariant ellipsoid, it is necessary to obtain a LDI of the saturation term as shown in the previous section (see Property 1). Therefore, the inequality (41) must be satisfied for the extremes of the convex hull, Ku¯x+du Lsand Hs¯x+h. 3.2.1. Case Ku¯x+du Ls Using Property 2 (40) and (28) in (39) is obtained that −vTPsv+2vTPs(ABc¯x+BBc Ku¯x+du Ls +EBc)− −λ1¯xTPs¯x+ρ(λ1−1) ≤0,∀¯x, ∀v The matrix representation of the previous inequality is given by νT⎡ ⎢ ⎢ ⎢ ⎣ ρ(λ1−1) ∗∗ 0−λ1Ps∗ −Ps(EBc+BLBcdu)−Ps(ABc+BLBcKu)−Ps ⎤ ⎥ ⎥ ⎥ ⎦ ν≤0(42) ∀¯x, ∀v, where BLBc =1 LsBBc,ν=1¯x−vT,and∗ represents the transpose of the symmetric term. So, ⎡ ⎢ ⎢ ⎢ ⎣ ρ(1 −λ1)∗∗ 0λ1Ps∗ Ps(EBc +BLBcdu)Ps(ABc +BLBcKu)Ps ⎤ ⎥ ⎥ ⎥ ⎦ >0 (43) Pre- and post-multiplying by diag[IP −1 sP−1 s] and making W=P−1 s ⎡ ⎢ ⎢ ⎢ ⎣ ρ(1 −λ1)∗∗ 0λ1W∗ (EBc +BLBcdu)(ABc +BLBcKu)WW ⎤ ⎥ ⎥ ⎥ ⎦ >0 (44) Considering that EBc =E+Axe+Brre−xe,BLBc = 1 LsBBc =1 LsBuand du=Cuxe+Dure−uc, it results that EBc +BLBcdu=Axexe+Brere+Ee where Axe=A−I+( 1 Ls)BuCu,Bre=Br+( 1 Ls)BuDuand Ee=E−(1 Ls)Buuc. On the other hand, reminding that ABc =A+BrKsand Ku=Cu+DuKs, it is obtained that ABcW+BLBcKuW=AwW+ByY where Aw=A+1 LsBuCu,By=Br+1 LsBuDu,andY= KsW. Therefore, the final LMI results as follows ⎡ ⎢ ⎢ ⎢ ⎣ ρ(1 −λ1)∗∗ 0λ1W∗ Axexe+Brere+EeAwW+ByYW ⎤ ⎥ ⎥ ⎥ ⎦ >0 (45) 3.2.2. Case Hs¯x+h In order to solve the case for Hs¯x+h,itiseasytosee that the same previous LMI is obtained only substituting Kuby Hs,duby h,andBLBc by BBc ⎡ ⎢ ⎢ ⎢ ⎣ ρ(1 −λ1)∗∗ 0λ1W∗ (EBc +BBch)(ABc +BBcHs)WW ⎤ ⎥ ⎥ ⎥ ⎦ >0 (46) Then, knowing that EBc =E+Axe+Brre−xe,BBc =Bu and ABc =A+BrKs: EBc +BBch=Anlxexe+Bnlrere+Bhh+Enle (ABc +BBcHs)W=AnlW W+BnlY Y+BvV where Anlxe =(A−I), Bnlre =Br,Bh=Bu,Enle =E, AnlW =A,BnlY =Br,Bv=Bu,Y=KsW,andV= HsW. So, the resulting LMI is ⎡ ⎢ ⎢ ⎢ ⎣ ρ(1−λ1)∗∗ 0λ1W∗ Anlxexe+Bnlrere+Bhh+Enle AnlW W+BnlY Y+BvVW ⎤ ⎥ ⎥ ⎥ ⎦ >0(47) For a fixed value of λ1, notice that the obtained matrix inequalities are LMI in the decision variables xe,re,ρ,Ps, Ks,Hs,andh. The procedure to find correct values for λ1 will be described later. Property 3 Suppose that there exists λ1≥0and λ2≥0 such that the LMI (35), (36), (45) and (47) are fulfilled. Then, ℘(Ps,ρ)is an invariant ellipsoid using the control law r=re+Ks(x−xe)and containing the system initial conditions, where P=W−1and K=W−1Y. Remark 1 Notice that the previous LMI depends multilineally on the parametric vector ˜ φ∈Φdue to the dependence of Ap=Ap(˜ φ),Bp=Bp(˜ φ),andEp=Ep(˜ φ). Then, properties 1 and 3 must be satisfied for all extreme plants of the hyperrectangle Φ. 3.3. Including performance inequality. Consider the representation of system (27) for the instant time k ¯x+ k=ABc¯xk+BBcσKu¯xk+du Ls+EBc (48) yk=Cy¯xk+Cyxe and suppose the following equality w=Cyxe(49) For an initial condition x0and the reference w, it is desired to calculate the system input rkby the law rk= 6 re+Ks(x−xe) such that the the following functional is minimized J= N  k=0 (yk−w)Q(yk−w)+¯xTK uRuKu¯x(50) where Q>0andRu>0 are symmetric matrices positive semi-defined. From the equality (49), it results that the functional J can be rewritten as J= N  k=0 ¯x kC yQCy¯xk+¯xTK uRuKu¯x(51) Defining LJ(¯xk)=¯x kC yQCy¯xk+¯xTK uRuKu¯x, it results that J= N  k=0 LJ(¯xk) (52) In the following property a strategy is defined for a correct selection of Ks,xeand rein order to fulfill the performance criteria (50). Property 4 Suppose that ¯x k+1Ps¯xk+1 −¯x kPs¯xk≤−LJ(¯xk)+γ, ∀˜ φ∈Φ,∀¯x and that an initial condition is equal to x0. Suppose also that the control law rk=re+Ks(xk−xe)is applied to the system, then J≤¯x 0Ps¯x0+Nγ where ¯x0=x0−xe. Proof 1 The assumption of the property leads to ¯x k+1Ps¯xk+1 −¯x kPs¯xk≤−LJ(¯xk)+γ, ∀˜ φ∈Φ,∀k≥0 Therefore, ¯x 1Ps¯x1−¯x 0Ps¯x0≤−LJ(¯x0)+γ ¯x 2Ps¯x2−¯x 1Ps¯x1≤−LJ(¯x1)+γ . . . ¯x NPs¯xN−¯x N−1Ps¯xN−1≤−LJ(¯xN−1)+γ ¯x N+1Ps¯xN+1 −¯x NPs¯xN≤−LJ(¯xN)+γ If the previous inequalities are added, it is obtained that ¯x N+1Ps¯xN+1 −¯x 0Ps¯x0≤−J+Nγ J≤¯x 0Ps¯x0+Nγ −¯x N+1Ps¯xN+1 J≤¯x 0Ps¯x0+Nγ So, from the previous property the following optimization problem can be proposed min Ps,Ks,xe,re,γ ¯x 0Ps¯x0+Nγ ∀˜ φ∈Φ,∀¯x(53) s.t. (¯x+)TPs(¯x+)−¯xPs¯x<−¯xC yQCy¯x−¯xTK uRuKu¯x+γ in order to calculate the control law that minimizes an upper limit of the functional. Notice that Nis a design parameter absent from problem constraints. If Nis very large, the problem solution will tend to minimize γ. The following interpretation can be considered for γ lim N→∞ J N≤lim N→∞ ¯x 0Ps¯x0+Nγ N=γ(54) Therefore, if Nis very large the initial transitory is almost not considered (the initial condition ¯x0doesn’t play a relevant role) and the emphasis is placed on improving the future behavior. On the other hand, if Ntakes very small values, the initial condition and the initial transitory gains relevance. Then, the problem (53) can be reformulated as min Ps,Ks,xe,re,γ,αs αs(55) s.a. ¯xT 0Ps¯x0+Nγ < αs (¯x+)TPs(¯x+)−¯xTPs¯x<−¯xTCT yQCy¯x− −¯xTKT uRuKu¯x+γ The problem inequalities will be translated to LMI form in order to address the optimization problem. Firstly, the upper inequality is considered ¯xT 0Ps¯x0+Nγ < αs(56) This can be easily expressed as a LMI using the Schur complement [25] in the form ⎡ ⎣ αs−Nγ ¯x(0)T ¯x(0) W⎤ ⎦≥0 (57) On the other hand, and remembering the presence of the saturation term in (48), the another inequality (¯x+)TPs(¯x+)−¯xTPs¯x<−¯xTCT yQCy¯x−¯xTKT uRuKu¯x+γ(58) must be satisfied for two extreme vertices of the LDI, Ks¯x+du Lsand Hs¯x+h, in the same way that for the invariant ellipsoid. The first step consists in using Property 2 on the previous inequality, where it is obtained that −vTPsv+2vTPs(ABc¯x+BBcσ(Ku¯x+du L)+EBc)−¯xTPs¯x< −¯xTCT yQCy¯x−¯xTKT uRuKu¯x+γ Following a similar procedure that in section 3.2, this inequality results in the following LMIs for the two extreme vertices of the LDI (the detailed procedure can be found in Appendix A [22]) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ γ∗∗∗∗ 0W∗∗∗ Axexe+Brere+EeAwW+ByYW∗∗ 0Q1/2CyW0I∗ 0RWW+RYY00I ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ >0 (59) 7 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ γ∗∗∗∗ 0W∗∗∗ Anlxexe+Bnlrere+Bhh+Enle AnlW W+BnlY Y+BvVW∗∗ 0Q1/2CyW0I∗ 0RvV00I ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ >0 (60) where Axe=A−I+(1 Ls)BuCu,Bre=Br+(1 Ls)BuDu, Ee=E−(1 Ls)Buuc,Aw=A+1 LsBuCu,By=Br+ 1 LsBuDu,Y=KsWand R1/2 uKuW=R1/2 u(Cu+DuKs)W= =R1/2 uCuW+R1/2 uDuKsW=RwW+RyY Finally, the equality (49), which was supposed before, must be included in the optimization problem. Therefore, the optimization problem has been reformulated to minimize the value of αssubject to a set of LMI. The following section describes the final optimization problem and the different obtained LMI. 3.4. Final optimization problem Notice that in section 2.0.2, it was considered that Ap= Ap(˜ φ), Bp=Bp(˜ φ), and Ep=Ep(˜ φ) for simplifying reasons. That is, it is necessary to remind that the matrices of the plant depend multilineally on the parametric vector ˜ φ. In this way, the previous LMI that were formulated for the nominal case, must be satisfied for all extreme values of the hyperrectangle Φ. Hence, the final problem can be formulated to calculate the decision variables xe,re,Ks, Ps,ρ,Hs,andh, in such a way that using the control law r=re+Ks(x−xe), ℘(Ps,ρ) is an invariant ellipsoid and the system fulfills the performance criteria given by J(50). The final optimization problem is given by min Ps,Ks,xe,re,γ,αs αs(61) s.a. ¯xT 0Ps¯x0+Nγ < αs (¯x+)TPs(¯x+)−¯xTPs¯x<−¯xTCT yQCy¯x− −¯xTKT uRuKu¯x+γ Then, considering the results obtained in previous section, a conservative way to solve the optimization problem consists in solving the following constraints w=Cyxe(62) ⎡ ⎣ αs−Nγ ¯x(0) ¯x(0) W⎤ ⎦>0 (63) ⎡ ⎣ ρ¯x(0) ¯x(0) W⎤ ⎦>0 (64) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ γ∗∗∗∗ 0W∗∗∗ Axe(˜ φ)xe+Bre(˜ φ)re+Ee(˜ φ)Aw(˜ φ)W+By(˜ φ)YW∗∗ 0Q1/2Cy(˜ φ)W0I∗ 0RW(˜ φ)W+RYY00I ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ >0(65) ⎡ ⎢ ⎢ ⎢ ⎣ ρ(1−λ1)∗∗ 0λ1W∗ Axe(˜ φ)xe+Bre(˜ φ)re+Ee(˜ φ)Aw(˜ φ)W+BYYW ⎤ ⎥ ⎥ ⎥ ⎦ >0(66) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ γ∗∗∗∗ 0W∗∗∗ Anlxe(˜ φ)xe+Bnlre(˜ φ)re+ +Bh(˜ φ)h+Enle(˜ φ) AnlW (˜ φ)W+BnlY Y+ +Bv(˜ φ)V W∗∗ 0Q1/2Cy(˜ φ)W0I∗ 0RvV00I ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ >0 (67) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ρ(1−λ1)∗∗ 0λ1W∗ Anlxe(˜ φ)xe+Bnlre(˜ φ)re+ +Bh(˜ φ)h+Enle(˜ φ) AnlW (˜ φ)W+BnlY Y+ +Bv(˜ φ)V W ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ >0(68) ⎡ ⎣ 1−h−λ2ρV V4λ2W⎤ ⎦≥0 (69) ⎡ ⎣ 1+h−λ2ρV V4λ2W⎤ ⎦≥0 (70) where it is necessary to incorporate constrains for each extreme value of the hypercube Φ. Also, as observed from the resulting constraints, some of them are BMI (Bilinear Matrix Inequalities) containing different bilinear terms ρ(1 −λ1), λ1W,λ2ρ), and 4λ2W. So, in order to obtain a stable MPC controller with good performance, it is necessary to choose λ1and λ2in a convenient way. Property 5 Suppose that there exist λ1≥0and λ2≥0 such that the constraints (62), (63), (64), (65), (66), (67), (68), (69), and (70) are feasible for every extreme of the 8 hypercube Φ. Then, there exits a control law r=re+Ks(x− xe)providing that ℘(Ps,ρ)is an invariant ellipsoid and the system fulfills the performance criteria given by J(50), where P=W−1and K=W−1Y. 3.5. Branch & Bound algorithm for bilinear terms Due to the presence of bilinear terms in the final optimization problem, a Branch & Bound algorithm is proposed to find the optimal solution. Before describing the algorithm, some aspects must be considered [22]: – From the constraints it is easy to see that λ1∈[λ1, λ1)= [0,1) and λ2∈[λ2, λ2)=[0,∞), respectively. Notice that λiand λirepresent the minimum and maximum values for λiwith i=1,2. –Alower bound solution is that obtained considering that there exist λ1∈[λ1, λ1)andλ2∈[λ2, λ2) such that the constraints (62), (63), (64), (65), and (67) are feasible, and also the following constrains are satisfied ⎡ ⎢ ⎢ ⎢ ⎣ ρ(1−λ1)∗∗ 0λ1W∗ Axe(˜ φ)xe+Bre(˜ φ)re+Ee(˜ φ)Aw(˜ φ)W+BYYW ⎤ ⎥ ⎥ ⎥ ⎦ >0(71) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ρ(1−λ1)∗∗ 0λ1W∗ Anlxe(˜ φ)xe+Bnlre(˜ φ)re+ +Bh(˜ φ)h+Enle(˜ φ) AnlW (˜ φ)W+BnlY Y+ +Bv(˜ φ)V W ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ >0(72) ⎡ ⎣ 1−h−λ2ρV VT4λ2W⎤ ⎦>0 (73) ⎡ ⎣ 1+h−λ2ρV VT4λ2W⎤ ⎦>0 (74) The list of lower bound solutions is called Ol. –Anupper bound solution is that obtained for specific values of λ1and λ2. Given a lower bound solution defined by [λ1, λ1) and [λ2, λ2), specific values for λ1and λ2are obtained as λ1=(λ1+λ1)/2andλ2=(λ2+λ2)/2. The list of upper bound solutions is called Pl. – The ranges [λ1, λ1) and [λ2, λ2) define an initial square which determines the search space. Such square will successively be divided until the optimal solution is obtained. Therefore, the proposed Branch & Bound algorithm presents the following steps: (i) Initialize Ol={},Pl={},[λ1, λ1)=[0,1), and [λ2, λ2)=[0,∞). (ii) Obtain the initial node noas an optimistic solution using the initial values for [λ1, λ1) and [λ2, λ2); Ol= {no}. (iii) The best node from the optimistic list Olis selected obtaining [λ1, λ1 ) and [λ2, λ2 ): (a) For odd iterations, the node will be that containing the square with bigger area. (b) For even iterations, the node will be that with better performance. (iv) New upper bound solution is calculated using the best node calculated in the previous step, [λ1, λ1 )and [λ2, λ2 ). The solution is calculated trying to solve the optimization problem for λ1=λ 1m=(λ1+ λ1 )/2andλ2=λ 2m=(λ2+λ2 )/2. If the optimization problem has solution for these values of λ1and λ2, the new value is included in the list of pessimistic solutions Pl. Then a local search is performed trying to improve the found solution. The local search is based on the following heuristic steps: (a) The values obtained for λ1and λ2, and the obtained performance solution are considered as initial variables. (b) Solve the optimization problem (61) using λ1 and λ2and obtaining values for Wand ρ.If a better performance solution is obtained, then the new upper bound solution is included in the upper bound list Pl. The best performance solution is updated with the new one. (c) Reformulate the LMI considering λ1and λ2as decision variables and, Wand ρas constant values. Solve the optimization problem obtaining new values for λ1and λ2. (d) Return to step (b) while better performance solutions are obtained or until a finite number of iterations is reached. (v) Remove worst lower bound solutions. The lower bound solutions with worst performance than the best upper bound one are removed from the lower bound list Ol. (vi) New lower bound solutions are obtained. That is, four new nodes are calculated from the node obtained in the previous step as follows n1={[λ1,λ  1m],[λ2,λ  2m]}, n2={[λ1,λ  1m],[λ 2m, λ2 ]}, n3={[λ 1m, λ1 ],[λ2,λ  2m]}, n4={[λ 1m, λ1 ],[λ 2m, λ2 ]}. The lower bound solution calculated for ni,i= 1, .., 4 will be removed if the solution is empty (the optimization problem is not feasible), the obtained performance is worst than the best upper bound solution, or the square defined by the associated ranges on λ1and λ2is too small. Otherwise, it will be included in the lower bound list Ol. 9