or Bilinear Matrix Inequalities*

European Journal of Control (2007)5:451–467 # 2007 EUCA DOI:10.3166/EJC.13.451–467

GPC Robust Design Using Linear and/or Bilinear Matrix Inequalities1 Jose Vincente Salcedo and Miranzo Martı´ nez Department of Systems Engineering and Control, Polytechnic University of Valencia, Camino de Vera, 14, Valencia, Spain

This article proposes an integrated methodology for the robust design of a generalized predictive controller (GPC) based on linear and/or bilinear matrix inequalities (LMIs/ BMIs). Firstly, a GPC formulation in state space which incorporates a full rank observer to estimate the on-line states is presented. It is then shown that the poles of such an observer are equivalent to the roots of filter polynomials Tj ðz1 Þ in a GPC Input/Output (I/O)formulation, which implies that these poles inherit the robustness effect associated with filter polynomial roots. In this way, we obtain closed-loop robust stability conditions and closedloop robust 1 -norm bounding conditions for the design of a GPC, based on LMIs and/or BMIs. These conditions take into account linear fractional (LFR) plant uncertainty. This type of uncertainty also contains affine dependence as a particular case, and can represent any matrix rational function. Finally, the robust GPC parameters are optimally obtained using a multiobjective genetic algorithm (GA). This integrated methodology is applied, as an example, to the multiobjective maximization of an uncertainty range and 1 -norm minimization in a benchmark plant: the two mass-spring mechanical system. Keywords: BMIs, GA, LMIs, predictive control, robust control.

1. Introduction The generalised predictive controller (GPC) was introduced by Clarke [7] as a particular case of model based predictive controllers. A CARIMA2 single input, single Correspondence to: J.V. Salcedo, Tel: þ 34-96387700 ext: 75766. Fax: þ 34-963879579. E-mail: [email protected] E-mail: [email protected] 1 Partially supported by the FEDER-CICYT research projects DPI2004-08383-C03-02 and DPI2005-07835 Spain.

output (SISO) I/O model is used for the process, and disturbances are assumed to be coloured noise. This methodology presents some advantages with respect to previous methodologies, as it enables work with unstable plants, an impossible task with methods based on the open loop step of a plant, or the impulse response [8]. In [26], the authors have developed a formulation for this controller in state space for multivariable plants. This formulation requires the use of an observer to estimate model states. Some works do not require the use of an observer, [5], although in many cases, a Kalman filter has been used [6]. Alternatively, in [35] a full rank observer using pole placement has been designed. After comparing the different alternatives, the use of a Kalman filter to estimate the states appears to be the most adequate choice. However, reference [24] states that it is difficult to select the variance-covariance matrices of the Kalman filter for precise tuning, specifically when using closed-loop specifications. For the reasons stated above, a full rank observer, designed by pole placement, is proposed by the authors in [26]. The design of the observer will take into account that the poles coincide with the roots of filter polynomials Tj ðz1 Þ which will be used when the GPC is formulated in the I/O version [5, 7]. This result implies that the observer poles inherit the robustness related to the filter polynomials Tj ðz1 Þ [22, 36]. Consequently, the design of this full rank observer is directly related to the closed-loop robustness, and as such, the poles of this observer should be selected accordingly when producing a robust GPC. In this article, the selection of observer poles and the remaining controller parameters, is made using linear 2 Controlled AutoRegressive Integral Moving Average. See [20]. Received 13 December 2006; Accepted 4 July 2007 Recommended by E. Camacho, D. Clarke

452

J.V. Salcedo and M.Martı´nez

x being the states of deterministic process; x the states related to noise inputs (); A, B and C a realization of the deterministic process; and , F1 and  are matrices related to the stochastic part of CARIMA model, known as noise matrices. This model is a specific realization in state space of the CARIMA I/O model used in [5]:

matrix inequalities (LMIs) [4], and bilinear matrix inequalities (BMIs) [10]. These matrix inequalities can take into account several specifications [28]: robust stability; robust satisfaction of norm bounds (1, 2 and 1); robust verification of time domain specifications (settling time, overshoot, etc.); and robust constraint satisfaction. These general conditions are adapted to the particular GPC formulation [26] used in this work, obtaining LMIs and/or BMIs which are numerically tractable [4,19]3. The selection of observer poles and the remaining controller parameters is performed optimally using a multiobjective GA, under the condition of verifying the LMIs, and/or BMIs, that assure the desired closed-loop robust performance achieved by the GPC. The steps presented above define the integrated methodology proposed in this work, which is the main contribution presented by the authors. The remainder of the article is structured as follows: in Section 2, the design of the GPC formulated in state space is briefly presented. Section 3 develops the LMIs and/or BMIs conditions, which guarantee the aforementioned robust specifications for the GPC. Section 4, presents the integrated robust design methodology for the GPC based on LMIs and/or BMIs and the multiobjective GA. Section 5 shows the control of a benchmark mechanical system using the robust GPC designed in state space. Finally, conclusions are presented in Section 6. Notation: Co(  ) represents the convex hull of a set. cols(  ) represents the number of columns of a matrix. Ker(  ) represents the kernel of a matrix. rows(  ) represents the number of rows in a matrix. size(  ) represents the size of a matrix. Span(  ) represents the vectorial space spanned by the columns of a matrix. z represents the forward operator in discrete time.

0 B B B B @

A1 ðz1 Þ 0 .. .



0

A2 ðz1 Þ    .. 0 .

0 .. .

0

10

y1

1

CB y C CB 2 C CB . C CB . C A@ . A

yn 0 0    An ðz1 Þ 1 1 1 1 10 B11 ðz Þ B12 ðz Þ    B1m ðz Þ u1 B B ðz1 Þ B ðz1 Þ    B ðz1 Þ CB u C 22 2m B 21 CB 2 C CB . C ¼B .. .. .. .. B CB . C @ A@ . A . . . . 0

Bn1 ðz1 Þ 0 1 T1 ðz Þ B 

B B þB B B @

0 .. . 0

 0 1

T2 ðz Þ 



0

 .. .

0



   Bnm ðz1 Þ 1 0 1 0  C 1 CB  C 0 CB 2 C CB . C B.C .. C @ A . C A . 1 Tn ðz Þ n

um



ð2Þ In a state space control technique, it is necessary to use an observer. Some authors [6], have proposed the use of a Kalman filter as an observer, due to its straightforward design. However, reference [24] affirms that the main problem with this observer is the choice of covariance matrices, and, as such, it is very constructive to develop observers which can be designed using a criteria closer to typical closed-loop specifications (settling time, overshoot). Other authors have used a full rank observer designed by pole placement [35]. Accordingly, an alternative full rank observer to estimate CARIMA model states (1) is proposed:

2. GPC Reformulation in State Space

x  ^ðk þ 1Þ ¼ ½A  FC ^ðkÞ þ BuðkÞ x þ FyðkÞ

To design the GPC in state space, the CARIMA model is used in [26]: 

xðkÞ þ BuðkÞ  xðk þ 1Þ ¼ A þ ðkÞ  yðkÞ ¼ C xðkÞ þ ðkÞ

3

 xðkÞ xðkÞ ¼   x ðkÞ   A  B A ¼ ; B ¼ 0 I 0   F1 ; C ¼ ½ C   ¼ I

ð3Þ

In a general case, the matrix inequalities obtained for a specific application are polynomial (PMI), much more complex than LMIs and BMIs [4].

ð1Þ

This full rank observer corrects the estimation of the state using the difference between the output xðkÞ, and the measurement of the outprediction, C^ put, yðkÞ.

453

GPC Robust Design with LMIs and/or BMIs

As previously established, the observer design consists of selecting the poles to ensure that the observer is stable and that its dynamics are faster than those of a closed-loop (controller þ process), in order to guarantee the convergence of the observer. For this observer, the dynamics are imposed by matrix  which has the following structure:  ¼ A  FC,      A  F1 ðC Þ ¼  F2 0 I ð4Þ   A  F1 C   F 1 ¼ F2 C I  F2

u



   u ðzÞ y ðzÞ ¼ MðzÞ ; yðzÞ uðzÞ 10 0 1 0 1 xðk þ 1Þ xðkÞ A B B CB B C B C @ y ðkÞ A ¼ @ C D D,p A@ u ðkÞ A C De, D yðkÞ uðkÞ

This conclusion may be of significance since it is known that filter polynomials, in particular their roots, give robustness to the GPC controller [22, 36]. With the use of the new state space CARIMA model (1) and the proposed full rank observer (3), the poles of this observer inherit this robustness. Later, in Proposition 4.1, the relation between robustness and observer poles will be rigorously established. These results are not surprising, since in [7], the filter polynomial is called the observer polynomial, due to its involvement in the prediction of future values4.

ð5Þ The uncertainty matrix belongs to the set:   ¼  ¼ diagð1 Ir1 , ...,r Irnr ,1 Ic1 , ...,nc Icnc ,1 , ...,nF Þ :  nr nc nF P P P i 2 R, j 2 C, Dl 2 CFl Fl , ri þ cj þ Fl ¼ n i¼1

The robust design of controllers can be established as the fulfilment of some specifications by the

yðzÞ ¼ Nðz,ÞuðzÞ;

4

¼

Although this fact was not justified in this reference.

A C

B D



 þ

l¼1

Block diagonal matrix  has three kinds of elements: i 2 R represents physical parameters of the plant having uncertainty, j 2 C and l 2 CFl Fl represent dynamical uncertainty, the first being structured and the second unstructured. The upper LFR interconnection of Fig. 1 can be simplified by eliminating u and y from (5), giving an uncertain dynamical relation between y and u with state space representation:

3.1. Plant Uncertainty



j¼1

ð6Þ

3. Robust Design Conditions Based on LMIs and/or BMIs



y

M(z)

closed-loop when challenged by the presence of disturbances, noise, and/or plant uncertainty. In this work, the uncertainty structure of the plant has been considered to be linear fractional (LFR) [9,10,29]. This type of uncertainty is obtained using an upper LFR interconnection [27] between a linear system M(z) and the uncertainty matrix  (Fig. 1). The linear system M(z) has the input/output and the state space representations:

Proposition 2.1: If ,  and F1 matrices, called noise matrices in the CARIMA state space model (1), are designed to establish the observer poles, then the roots of the filter polynomials (Tj ðz1 Þ), corresponding to the equivalent input/output (I/O) version of the GPC controller are equal to the assigned observer poles.

xðk þ 1Þ yðkÞ

yΔ

Fig. 1. Upper linear fractional (LFR) interconnection providing a LFR uncertainty.

As a consequence, it is possible to select F1, F2,  and  to set the eigenvalues of  at the desired placements, thereby obtaining the corresponding dynamics of the observer. In particular, it is shown that F2 can be taken as the identity matrix without any loss of generality. An important result of [26] is:



⌬

uΔ

  B ðI  D Þ1 C D e, 

D,p





xðkÞ uðkÞ

 ð7Þ

454

J.V. Salcedo and M.Martı´nez

Given any fixed value for i, j and l this expression can transform into a linear system. In particular, the dependence with respect to these parameters is rational due to the matrix inversion. This LFR uncertainty is quite general and can be applied to real cases [4,10,17,23,31]. Moreover, the LFR uncertainty has, as a particular case, affine uncertainty – as can be seen in (7) taking D ¼ 0. 3.2. Linear and Bilinear Matrix Inequalities LMIs are very useful in robust control theory [4, 10, 27]. Specifically, they can be used for analyzing closedloop properties for a given controller, or for designing a robust controller to meet some robust specifications. An LMI is defined as [4]: ~

FðxÞ ¼ F0 þ

m X

xi Fi > 0,

5

ð8Þ

i¼1

Where x 2 Rm is the unknown vector and the symmetric matrices Fi ¼ FTi 2 Rnn , i ¼ 0, . . . , m are given.

Robust stability is the most commonly used specification in robust control. It is based on the Lyapunov stability condition which, for discrete linear time invariant (LTI) systems, can be formulated as: A discrete LTI system is asymptotically stable if, and only if, 9P ¼ PT > 0 in such a way that: ð9Þ

This condition is also called quadratic stability [4]. This inequality can be viewed as an LMI if the nonrepeated elements of the P matrix are taken as unknowns. When the plant has uncertainty, condition (9) cannot be applied. Instead, under uncertainty, the following quadratic stability condition can be used: 9P ¼ PT > 0 : AðÞT PAðÞ  P < 0

8 2 1 ð10Þ

This quadratic stability condition is only sufficient [4] for guaranteeing asymptotic stability under uncertainty. The set 1 represent all the possible values of the The inequality symbol > means that F(x) is a positive definite matrix. By definition, the previous LMI is strict, although it is possible to consider nonstrict LMIs using  instead of > . If the set defined by an LMI is nonempty then it is convex.

5

^ ¼ PðÞ ^ T>0: 9PðÞ ^ ^ < 0 8 2 1,8 ^ 21 ^ AðÞT PðÞAðÞ  PðÞ ð11Þ Henceforth, this dependence will be assumed to be LFR since it represents a general case of dependence. ^ Note that in the P matrix, the uncertainty matrix  can differ from  [16]. In actual fact, both can be dependent upon the same parameters, but probably with different sizes, for example:  ¼ diagð1 Ir1 , . . . ,r Irnr ,1 Ic1 , . . . ,nc Icnc ,1 , . . . ,nF Þ

ð12Þ

3.2.1. Robust stability

AT PA  P < 0

uncertainty matrix . The main problem of condition (10) is that it represents an infinite number of LMIs, each one for a fixed value of matrix . To overcome the problem of a infinite number of LMIs in (10) some methodologies have been developed in order to acquire an equivalent condition based in a finite number [13,29,30]. It is possible to acquire a less conservative condition for robust asymptotic stability using an uncertain dependent Lyapunov matrix [13, 16] in the quadratic stability condition with uncertainty (10):

 ^ ¼ diag 1 Ir0 , . . . ,r Ir0 ,1 Ic0 , . . . ,nc Ic0 ,1 , . . . ,nF  nr nc 1 1

ð13Þ where any of sizes rj, rj0 , cj and cj0 can be read as zero. For a plant with LFR uncertainty, and with LFR dependence on the Lyapunov matrix, it is possible to derive a finite number of LMIs which provide a sufficient condition for quadratic stability (11). As for a plant with LFR uncertainty: ~

AðÞ ¼ A þ B ðI  D Þ1 C

ð14Þ

A, B, D and C being constant matrices of proper dimensions. The Lyapunov matrix has an LFR dependence [16]:

1 ~ ^ I  TD  ^ ^ ¼ TA þ TB  TC TðÞ ð15Þ ~ ^ ¼ TðÞ ^ T QTðÞ ^ > 0 8 ^ 21 ^ PðÞ ^ is fixed a priori following the The structure of TðÞ ideas of [16]:  Taking TB ¼ 0, TC ¼ 0 and TD ¼ 0, an uncertain independent Lyapunov matrix (constant) is recov^ ¼ TT QTA . ered: TðÞ A    0 R S , TB ¼ , partied con Taking Q ¼ I ST 0 ^ has an affine depenformably and TD ¼ 0, TðÞ ^ dence on .

455

GPC Robust Design with LMIs and/or BMIs

^ has a quadratic dependence  Taking TD ¼ 0, TðÞ ^ on .  If all the matrices, TA ¼ 0, TB ¼ 0, TC ¼ 0 and ^ has a TD ¼ 0, are not equal to zero matrix, TðÞ ^ general LFR dependence on . ^ the only unknown is the Q Consequently, in PðÞ ^ matrix. PðÞ has, in general, an LFR dependence ^ and TðÞ ^ has, in gensince Q does not depend on  ^ In eral, an LFR dependence with respect to . Appendix A.2 a sufficient condition for (11) is developed: there is a symmetric matrix V such that: ! Q 0 0 T T T<0 ð16Þ 0 Q 0 V ~  IsizeðÞ ~

!T V

~  IsizeðÞ ~

! >0

ð17Þ

~ 2 X ~ , V11 < 0 8  where V11 is the (1,1)-block of symmetric matrix !V ~  . partied conformably integrated with IsizeðÞ ^ Finally, the robust stability problem can be recast as r þ 2 LMIs given by (16) and (17), where the unknowns are Q and V, since the uncertainty set X ~ is a priori knowledge6. The matrix T is obtained in Appendix A.2. Remark 3.1: These conditions are obtained in the case ~ is given as the convex hull of ~ (where belongs ) when 1 the finite set of matrices X ~ (see Remark A.1). Suffi~ are cient conditions under another hypothesis for 1 provided in Remarks A.2 and A.3.

LMIs (16) and (17) are presented. These techniques could be used for concluding when these sufficient LMIs are nonconservative. This paper also presents other relaxation techniques for obtaining a finite number of LMIs: the Po´lya relaxation, and sum of square decompositions. These methods are guaranteed, under some conditions, to be asymptotically exact. One disadvantage of these methods is that the degree of relaxation which gives an exact result has no known a priori bound. Such techniques could be used in cases when LMIs (16) and (17) do not provide a feasible solution.

3.2.2. The determination of a maximum uncertainty set which guarantees robust stability Robust stability conditions (16) and (17) can also be used to determine the maximum range of some real parametric uncertainties which guarantee robust stability, instead of checking robust stability for a priori given range. In this case, some matrices of set X ~ will have some entries as unknowns, hence the corresponding matri~ in the LMIs (50). Consequently, these matrix ces  inequalities become non-LMIs – due to products of up to three unknowns8: 

7

See Appendix A.2. Depending on particular uncertainty structure (6).

V11 VT12

V12 V22



~  I

 >0

~ 2X~ 8 

By incorporating LMI (51) and conducting some operations: V11 > 0, 1 ~T ~ ~ ~ T V12 þ VT   12 þ V22   V11 ðV11 Þ V11  > 0 ~ 2X ~ 8 

ð18Þ and applying the Schur complement (appendix A.1) at this point: ! ~ þ V22  ~ T V11 ~ T V12 þ VT   12 >0 ~ V11  V11 ~ 2 X~ 8 

Remark 3.4: In the survey paper [30], techniques for verifying the exactness approach of the conservative 6

T 

~ 2X~ ~ T V11  ~ þ ~ T V12 þ VT  ~ þ V22 > 0 8  12 

Remark 3.2: Instead of using the Lyapunov paradigm for addressing robust stability, standard -theory [23,37] could also used. In such a case, the upper bound computations for the complex, mixed, or real structured singular value7 must be performed over a frequency grid, instead of solving the set of LMIs (16) and (17). Remark 3.3: Also, the robust techniques of Ben-Tal and Nemirovski [2,3] could be used for obtaining a finite number of sufficient LMIs for condition (11), different from conditions (16) and (17), which are based on the ideas of work [29].

~  I

ð19Þ ~ Elements of V and .

8

456

J.V. Salcedo and M.Martı´nez

these non-LMIs become BMIs. A BMI [34] is a generalization of an LMI (8), incorporating products between two unknowns: m X

m X m X

In Appendix A.3 the LMI conditions (59) for 1norm bounding are deducted.

ð20Þ

3.2.4. Maximizing the uncertainty set and minimizing the 1-norm bound

The set defined by a BMI can be nonconvex, and the related problems related can be NP-hard [32]. In (19), ~ ~ and V11, and between  there are products between  and V12. Finally, the determination of the maximum uncertainty set can be recast as a optimization problem (the maximization of some uncertainty ranges) subject to LMI (16) and r BMIs (19), where the unknowns are Q, V and some entries for matrices belong to set X ~ .

It is possible to maximize the uncertainty set and minimize the 1-norm bound simultaneously using the matrix inequalities (59). This problem is a multiobjective optimization problem subject to (59). As seen in section 3.2.2 the matrix inequalities: !T ! ~ ~   ~ 2 X~ ð25Þ V > 0 8  IsizeðÞ IsizeðÞ ~ ~

~

FðxÞ ¼ F0 þ

xi Fi þ

i¼1

xi xj Fi,j > 0:

i¼1 j¼1

Remark 3.5: This idea could be used for other types of uncertainty, not only for real parametrics. For example, ~ has a full complex structure (see Remark A.2) if 1 ~ k : norm-bounded k1 1        ~ ¼  2 Cmn : T I I 0  0 1 0 2I I ð21Þ Consequently, the maximum size of this full complex uncertainty set is obtained maximizing  2 subject to (16) and conditions (see Remark A.2):   I 0 ð22Þ y11 0, V ¼ y11 0 2I Note that (16) becomes a BMI when substituting the multiplier V by the expression given by (22) due to the product between y11 and  2. Furthermore, the first condition of (22) is an LMI. 3.2.3.

1-norm bounding

Remark 3.6: Using the same ideas contained in Remark 3.5, it is possible to maximize the size of other types of uncertainty sets.

4. GPC Robust Design 4.1. Introduction The state space GPC control law is obtained solving the optimization problem: Min uðkÞ,,uðkþNu 1Þ

1-norm bounds related to an input/output channel can be expressed also in terms of a set of LMIs. In particular, to robustly bound the 1-norm of channel !/y by  for the linear uncertain plant: xðk þ 1Þ ¼ AðÞxðkÞ þ B! ðÞ!ðkÞ yðkÞ ¼ Cy ðÞxðkÞ þ Dy! ðÞ!ðkÞ

^ containing some are nonlinear since matrices  unknowns. In Section 3.2.2 these matrix inequalities were converted to BMIs (19). Consequently, when the objective is to maximize the uncertainty set and minimize the 1-norm bound, a multiobjective optimization problem must be solved subject to the first LMI of (59) and r BMIs (19), where the unknowns are Q, V,  and some entries for matrices belong to set X ~ .

ð23Þ

it is sufficient that the following infinite set of LMIs is satisfied by certain dependent Lyapunov matrices [12,25,27]:

^ ^ þ Cy ðÞT Cy ðÞ  PðÞ AðÞT PðÞAðÞ T ^ B! ðÞ PðÞAðÞ þ Dy! ðÞT Cy ðÞ

X N2

ð!ðk þ iÞ  yðk þ ijkÞÞT Si ð!ðk þ iÞ

N1 ¼1 Nu X

 yðk þ ijkÞÞ þ

 uðk þ iÞT Ri uðk þ iÞ

ð26Þ

j¼0

N1 and N2 being the beginning and the end of the prediction horizon, Nu the control horizon, y(k þ ijk) the output prediction at instant k þ i with the information available at instant k, !(k þ i) is the reference tracking, and Si and Ri are weighting matrices with proper dimensions.

^ ! ðÞ þ Cy ðÞT Dy! ðÞ AðÞT PðÞB T ^ ! ðÞ   2 I Dy! ðÞ Dy! ðÞ þ B! ðÞT PðÞB

! <0

^ 8 , ð24Þ

457

GPC Robust Design with LMIs and/or BMIs

From the optimal solution of this problem, only the control action corresponding to first instant in the control horizon, u(k), is applied to the plant. This process is repeated at instant k þ 1. This technique is known as receding horizon. In [26], a linear state space representation is obtained for this receding horizon control law which depends on observed states of CARIMA model (1): xðkÞ xu ðk þ 1Þ ¼ ðI  OÞxu ðkÞ  ðI  OÞM þ ðI  OÞ!ðkÞ  ðI  OÞPyðkÞ

Fig. 2. Block diagram for HðzÞ and RðzÞ generalized predictive controller (GPC) decomposition.

xðkÞ þ !ðkÞ  PyðkÞ uðkÞ ¼ Ixu ðkÞ  M ð27Þ and combining it with the proposed full rank observer (3), and so a dynamic output feedback controller is obtained: ~



xc ¼ x ðk þ 1Þ ¼



xu  A  C  BM

c

B

!

ðI  OÞM I  O |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} 

þ

x

 BU

Ac





ðI  OÞU |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}

!ðkÞ þ

Bc!

xc ðkÞ

   BP



ðI  OÞP |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

yðkÞ

Bcy

uðkÞ ¼ ð M I Þ xc ðkÞ þ |{z} U !ðkÞ |ffl{zffl} P yðkÞ |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} c c Cc

D!

Dy

ð28Þ , M, P, O and U are matrices which nonlinearly depend on weighting matrices, control and prediction horizons, and observer poles. Consequently, to design the GPC controller in state space it is necessary to select the following parameters:  Horizons: N1, N2 and Nu  Weighting matrices: Si and Ri  Observer poles. The objective of this section is to propose a methodology for selecting these parameters in order to robustly guarantee certain closed-loop specifications. In general, all the matrices of a GPC controller [26] Ac, Bc! , Bcy , Cc, Dc! and Dcy have a nonlinear dependence with respect to GPC design parameters.

4.2. Nominal Robust Analysis First of all, a robust analysis is carried out using sensitivity functions when plant uncertainty is not present. As such, only disturbances are considered. GPC can be viewed as a two degrees of freedom controller, as shown in Fig. 2 shows, where the HðzÞ

transfer matrix represents the part related to exogenous reference and the RðzÞ transfer matrix represents the part related  to output  measurement. In the same A B represents the nominal figure G0 ðzÞ ¼ C 0 plant, du and dy represent additive input and output disturbances, and dm noise measurement. By making some transformations in Fig. 2, a new equivalent block diagram is obtained, which has the structure of a classic multivariable closed-loop with controller RðzÞ and prefiltering transfer matrix R1 ðzÞHðzÞ. Taking the nomenclature of [15], the closedloop can be described using the transfer matrices: S0(z) nominal sensitivity, T0(z) nominal complementary sensitivity, and Si0(z) nominal input disturbance sensitivity. This provides the following relationship [15]: yðzÞ ¼ T0 ðzÞR1 ðzÞHðzÞ!ðzÞ  T0 ðzÞdm ðzÞ þ S0 ðzÞdy ðzÞ þ Si0 ðzÞdu ðzÞ ð29Þ These sensitivity functions can analyze some aspects of closed-loop performance [15]: reference tracking, perturbation rejection, noise measurement immunity, control effort, robustness against uncertainty and/or unmodelled dynamics, etc. Closed-loop transfer matrix can be obtained starting from (29): yðzÞ ¼ ðI  G0 ðzÞRðzÞÞ1 G0 ðzÞHðzÞ!ðzÞ

ð30Þ

Since in this analysis there is no uncertainty present in the plant, the separation principle is applicable and as such, it can be concluded that this transfer matrix is independent from the selection of observer poles, although it depends on horizon sizes, and weighting matrices. This result reveals that it is possible to select horizon sizes and weighting matrices to achieve some closed-loop specifications: settling time, bandwidth, etc., independently from observer pole selections.

458

J.V. Salcedo and M.Martı´nez

However, sensitivity functions depend on the selection of observer poles. The following proposition states this idea in a technical manner: Proposition 4.1: (T0(z) depends on observer poles)The nominal complementary sensitivity function, T0(z), has as poles closed-loop poles and observer poles – unless pole/zero cancellations are present. The proof of this proposition is outlined in Appendix A.5. Similar conditions apply for the remaining sensitivity functions. This result proves rigourously the influence of observer poles on closed-loop robustness. Therefore, a correct selection of observer poles can provide a closed-loop with some level of robustness. In the following subsection, a methodology based on LMIs and/or BMIs will be established to enable a selection of observer poles which robustly achieve some closedloop properties. 4.3. Integrated Robust Design Methodology Given the uncertain LTI discrete system subject to disturbances and noise measurement of Fig. 2, which represents the uncertain plant: xðk þ 1Þ ¼ AðÞxðkÞ þ BðÞðuðkÞ þ du ðkÞÞ yðkÞ ¼ CðDÞxðkÞ þ dy ðkÞ; ym ðkÞ ¼ yðkÞ þ dm ðkÞ ð31Þ

Once this closed-loop expression is available, to design a robust GPC with the LMI and/or BMIs conditions suggested in Section 3, for example (16) and (17) for robust stability, (16) and (19) for maximizing an uncertainty set, (59) for 1-norm bound in some input/output channel, or the first LMI of (59) and (19) for simultaneous maximization of uncertainty set and minimization of 1-norm, it is only necessary to find a set of GPC paramenters: horizons, weighting matrices and observer poles that are capable of satisfying these matrix inequalities (16), (17), (19), and (59). Due to the nonlinear dependence of closed-loop matrices (28) with respect to GPC parameters and the possible use of robust BMI conditions, the selection of GPC parameters becomes a nonconvex problem. Consequently, in this work, this set of parameters is obtained using a multiobjective genetic algorithm (GA) [21] since it is a global optimization method. In particular:  The GA searches through the set of all the possible combinations of GPC parameters, taking these as individuals of the population.  The GA sorts the individuals from the population that satisfy the LMIs and/or BMIs, using a fitness function. The fitness function is dependent upon the robust problem which must be solved. In general, the formulation of the GA problem can be stated as: min

½f1 , . . . ,fn 

GPCparameters

subject to the LMIs and/or BMIs

the closed-loop with GPC is (28):

ð33Þ

Acl

zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl3{ 2 AðÞ  BðÞPCðÞ BðÞM BðÞ     xðk þ 1Þ 6 7 xðkÞ      ¼4 ð  BPÞCðÞ A  C  BM B 5 c þ xc ðk þ 1Þ x ðkÞ ðI  OÞPCðÞ ðI  OÞM I  O Bcl,dym Bcl,! zfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflffl3 ffl{ zfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflffl{ Bcl,du 2 2 3 zfflfflfflfflffl BðÞP BðÞP  ffl}|fflfflfflfflfflffl{ BðÞ 6 7 6 7   du ðkÞ þ 4   BP þ4 BU 5 ðdy ðkÞ þ dm ðkÞÞ 5 !ðkÞ þ 0 ðI  OÞP ðI  OÞU Ccl,y D cl,yd zfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflffl{  xðkÞ  z}|{ yðkÞ ¼ ½ CðÞ 0  c þ I dy ðkÞ x ðkÞ Ccl,u Dcl,ud Dcl,u! zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{  xðkÞ  z}|{ zffl}|ffl{ þ U !ðkÞ P ðdy ðkÞ þ dm ðkÞÞ uðkÞ ¼ ½ PCðÞ M I  c x ðkÞ

ð32Þ

459

GPC Robust Design with LMIs and/or BMIs

Each fj depends on plant uncertainty, GPC parameters and robust specifications. As such, the fitness function is ½ f1, . . . , fn. The solution for multiobjective optimization problems will be an approximation of the optimal Pareto front composed by a finite number of individuals from the population (GPC parameters). For simplicity of notation, the matrix inequality which groups all the LMIs and/or BMIs of a given problem, is introduced: F(plant uncertainty, GPC parameters, robust specifications) < 0

ð34Þ

 Given an individual or GPC parameters, the GA first check if these satisfy the matrix inequality F < 0. If so, then they become valid GPC parameters and the GA evaluates their corresponding fitness value. In particular cases, the fitness function can be taken as: 1. If the objective is to maximize the set of allowable uncertainties, then the fitness function can be taken as a certain size measure of that set. For example, if the uncertain plant has an uncertain real parameter, , whose range of variation has to be maximized, then f1 ¼ ðmax  min Þ. This case is single objective. 2. If the objective is to minimize the 1-norm of an input/output channel (59) then the performance function will be this 1-norm value: f1 ¼ . 3. If the objective is to maximize an uncertainty set and minimize the 1-norm of an input/output channel, then, the fitness function will be: ½f1 ,f2  ¼ ½ðmax  min Þ,. The checking of matrix inequality F < 0 is performed using two methods. If all the inequalities are LMIs, there are very efficient algorithms which can solve these in polynomial time [4,10]. Some packages with Matlab can solve these, for example: LMI toolbox distributed by Mathworks [14] or SeDuMi9 which is free. However, in general, F > 0 will contain BMIs and, as such, the problem will be nonconvex and NP-hard [32,34]. For these reasons, there are no general numerical algorithms which can solve these in polynomial time. Some authors [11,33] have proposed global algorithms to solve them. These algorithms have a high computational cost and, in practice, are applicable only to small problems. Other authors [19], have proposed local algorithms that only reach local 9

Available at http://sedumi.mcmaster.ca.

optima, but which can handle large BMIs and, furthermore, the computational cost is low compared with global algorithms. There is commercial software that implements the algorithm proposed in the previous reference: PENBMI10. This software will be used here to locally solve the BMI problems that occur. A potential weakness of the algorithm described above is that the number of GPC parameters to be fitted and that the size of LMIs and/or BMIs grows with the complexity of the plant (number of states, number of inputs and number of outputs). This implies that the time required for solve the GA algorithm with a fixed population and a fixed number of generations could be very expensive. Instead, nonlinear programming techniques could be used as they require a less computation effort. However, the GA guarantees that the global optimum will be pursued without stalling at discontinuity points and avoiding local optimums, as opposed to nonlinear programming which reach local optimum only. All the computations are made off line, and consequently, the computing requirements could be relaxed.

5. Application Example This example is a process benchmark used in [18]. It is a mechanical system formed by two masses and one spring (Fig. 3): 2 3 0 12 3 0 0 1 0 x1 x_ 1 6 x_ 7 B 6 C 0 0 0 1 C6 x2 7 6 27 B 7 6 7¼B C6 7 4 x€1 5 @ K=m1 K=m2 0 0 A4 x_ 1 5 x€2

K=m2 1

0

K=m2

0 B 0 C B C þB Cu, y ¼ ð 0 @ 1=m1 A

0

1

0

0

x_ 2 3 x1 6x 7 6 27 0 Þ6 7 4 x_ 1 5

0

2

x_ 2 ð35Þ

where m1 ¼ m2 ¼ 1 kg and K 2 ½0.5 10 N/m, that is to say, the spring constant has uncertainty. The nominal value of K is taken as the mean value of uncertain range: 5.25. With this notion: KðÞ ¼ 5:25 þ 4:75  , jj 1

ð36Þ

This process sets its four poles at the imaginary axis, and the controlled variable, y ¼ x2, is not directly affected when using the control action, u. For these reasons [18], designing adequately robust stabilizing 10

www.penopt.com.

460

J.V. Salcedo and M.Martı´nez

0

x u

1 0 B 0 C C B~d ¼ B @ T=m1 A, Cd ¼ ð 0 0

y = x2

K

m

m

Fig. 3. Two mass-spring mechanical system

controllers is a demanding task. In particular, the process has affine uncertainty with respect to the uncertain spring constant K. The inclusion of more uncertain parameters, such as m1 or m2, provides a LFR uncertainty instead. However, this makes the problem even harder to solve. As such, in this paper only K will be considered as the uncertain parameter. First a discrete time model is obtained applying a zero order holder [1]: Z T Ac T , Bd ¼ eAc t Bc dt ð37Þ Ad ¼ e 0

T being the sampling period, Ac and Bc the continuous time model, and Ad Bd the discrete model. If the continuous model has uncertainty, the discrete one will too, but in a more complex way due to the matrix exponential. Specifically, if the continuous time model has affine uncertainty, the discrete model will have exponential uncertainty. In this study, the most commonly proposed uncertainty is LFR, which is less complex than exponential uncertainty. For this reason, an approximate method to compute (37) based on a series expansion, using the exponential series, for discrete matrices is proposed11. Taking a small sampling period compatible with hardware specifications and controller computational requirements, only a finite number of terms corresponding to this expansion will be used to guarantee a maximum error between exact and approximate discretizations. Since the finite series has a polynomial dependence with respect Ac, if Ac uncertainty is LFR then Ad uncertainty will also be LFR. In addiction, if Bc has LFR dependence (37) then Bd will also have LFR dependence. In [18] a sampling period of T ¼ 0.1 s and a series expansion of order one in T for both discrete matrices is used, or equivalently, the Euler discretization. In this study, a smaller sampling period T ¼ 0.01 and a first order expansion are also used: 0 1 1 0 T 0 B 0 1 0 TC B C A~d ¼ B C, @ KT=m1 KT=m1 1 0A KT=m2 11

KT=m2

Bilinear approximation can also be used.

0

1

1 0

0Þ

ð38Þ

This sampling period is adopted to achieve a maximum norm error between exact and approximate matrices of 0.001. For the nominal value of K, the norm errors between state matrices and input matrices are: keAc T  A~d k2 ¼ 5:3336  104 Z T k eAc t Bc dt  B~d k2 ¼ 5:0013  105

ð39Þ

0

and for the remaining values of K corresponding to interval ½0.5, 10, the maximum error of 0.001 is always satisfied. To reiterate, when a finite number of terms in the series expansion is used, and if the continuous time uncertainty is LFR, then the discrete model will also be LFR. As the continuous model for the two massspring process has affine uncertainty(35), and a first order series expansion in T is performed for both discrete matrices (38), then the approximate discrete time model also has affine uncertainty, a particular case of LFR. The next step is to establish the GPC design requirements:  Maximize range of variation of K, min ¼ max  max , around its nominal value with guaranteed closed-loop robust stability.  Minimize 1-norm, , corresponding to reference/y channel. In general, these conditions are contradictory and they provide a multiobjective nonconvex optimization problem, which will be solved using a multiobjective GA. The multiobjective GA employs the following ranges for GPC parameters (individuals from the population): 1 N2 256 1 Nu 8 Observer poles 2 ½0, 0.99 Error weighting matrices: Qi ¼ 1 Control action weighting matrices: 0 l 100  K() ¼ 5.25 þ 4.75  

    

and the following genetic parameters:  Population size: 2000  Crossover probability: 0.8

Ri ¼ l,

461

GPC Robust Design with LMIs and/or BMIs

20

12 11

δ = –0.83 0

10 9

–20 Gain (dB)

∞−Norm

8 7 6 5

–40

–60

4 3

–80

2

1.1

1 0

–100 100

0.8017

101 Frequency (rad/s)

0

0.2

0.4

δmax

0.6

0.8

1

Fig. 4. Approximate optimal Pareto front for the multiobjective optimization problem related to generalized predictive controller (GPC) robust design for a two mass-spring system.

 Mutating probability: 0.1  Multiobjective fitness function: equivalent to ½  max, .

½  2max,



With a pentium IV 2.8 GHz under Windows XP using Matlab R2006a, the average time for computing a generation was about 2 hours. This time can be efficiently reduced using the following ideas: programming the code in C language instead of using Matlab, and using a parallel distributed computation based on a multiprocessor system. In Fig. 4, an approximation to the optimal Pareto front corresponding to the multiobjective optimization problem obtained following 300 generations with GA is shown. This is formed by a finite number of individuals, each one corresponding to some GPC parameters. In some cases, there are gaps between these, which means that the GA did not obtain any individuals belonging to this zone of the Pareto front. It would be necessary to run many more generations using the GA in an attempt to complete these gaps. In this figure, it can be seen that for lower values of max, 1-norm is approximately 1, while the maximum value of max guaranteeing closed-loop stability is 0.916 with a 1-norm of 10.9001. Both situations result in an extreme GPC controller: the first with a minimum 1-norm, and the second with a maximum range of K providing robust closed-loop stability. From this figure, GPC controllers can be obtained which verify:  GPC with minimum 1-norm for a given maximum uncertainty range, max.

Fig. 5. Closed-loop Bode diagrams for  0.83  0.9.

 GPC with maximum uncertainty range, max, for a given 1-norm upper bound. For example, if the GPC controller with a maximum uncertainty range for an 1-norm bounded by 1.1 was obtained, using the Pareto front (Fig. 4) it could be deduced that max ¼ 0.8017 and the GPC parameters would be:  Nu ¼ 3, N2 ¼ 164, Ri ¼ l ¼ 0.0046.  Observer poles: 0.2049, 0.2800, 0.0530, 0.0368 and 0.3339. Their equivalent filter polynomial is Tðz1 Þ ¼ 1  0:9086z1 þ 0:29476z2  0:040445z3 þ 0:0021479z4  3:7363  105 z5 : this max implies that maximum range of variation of K is ½1.4420, 9.0580. By computing the GPC corresponding to these parameters, and simulating the closed-loop for different values of , it can be concluded that the closed-loop is robustly stable, and approximately in the range  0.83  0.9, that is to say, 1.3075 K 9.525. However, in this range there is no guarantee that the 1-norm will be bounded by 1.1, specifically, for K ¼ 1.3075 ( ¼  0.83), the 1-norm is 1.3411. In Fig. 5, closed-loop Bode diagrams for values of  belonging to interval ½  0.83, 0.9 are represented. The largest Bode diagram peak is obtained when  ¼  0.83. In Fig. 6 unit step closed-loop responses corresponding to  ¼ 0, 0.9 and  0.83 are represented, and it can be seen that case  ¼  0.83 has a worse performance than other cases (longer settling time), while the others have a very similar response.

462

J.V. Salcedo and M.Martı´nez

References 0.7

0.6

Output & reference

0.5

0.4 δ = −0.83 0.3

0.2

0.1

0

0

200

400

600

800

1000

Samples

Fig. 6. Unit step closed-loop responses for cases  ¼ 0, 0.9 and  0.83.

6. Conclusions 1. This paper presents an integrated robust design methodology for the GPC. It consists of the optimal selection of GPC parameters using a multiobjective GA, under the condition of verifying the LMIs and/or BMIs that assure the desired closedloop robust performance achieved by the GPC. 2. To achieve this, it is initially and formally established that the selection of observer poles can provide some robust properties for the closed-loop. The proof is conducted using the sensitivity functions. 3. Then robust stability and 1-norm conditions based on LMIs and/or BMIs for the GPC formulation presented in state space are obtained. 4. Finally, a multiobjective GA algorithm is used to robustly design the GPC. The fitness function of this GA depends on the robust specifications to be achieved. Inside the GA, the numerical methods to solve LMIs and BMIs are utilized. 5. The integrated methodology for designing robust GPCs is applied to a benchmark process: the two mass-spring system. With this design, a multiobjective optimization problem is solved: maximizing the uncertainty size of K versus minimizing the 1-norm of the reference/output channel.

Acknowledgements We would like to thank the R þ D þ i Linguistic Assistance Office at the Universidad Politecnica of Valencia for their help in revising this paper.

1. Aström KJ, Wittenmark B. Computer controlled systems. Theory and Design. Prentice Hall Information and Systems Sciences Series, 1997 2. Ben-Tal A, ElGhaoui L, Nemirovski A. Robust Semidefinite Programming. Handbook on Semidefinite Programming, Kluwer Academic Publishers, 2000; pp. 139–162 3. Ben-Tal A, Nemirovski A. On tractable approximations of uncertain linear matrix inequalities affected by interval uncertainty. SIAM J Optim 2002; 12: 811–833 4. Boyd S, El Ghaoui L, Feron E, Balakrishnan V. Linear Matrix Inequalities in System and Control Theory. SIAM Studies in Applied Mathematics, Philadelphia, 1994 5. Camacho EF, Bordons C. Model Predictive Control. Springer, London, 2004 6. Chisci L, Lombardi A, Mosca E, Rossiter J. State-space approach to stabilizing stochastic predictive control. Int J Control 1996; 65(4): 619–637 7. Clarke D, Mohtadi C. Properties of generalized predictive control. Automatica 1989; 25(6): 859–875 8. Cutler C, Ramaker D. Dynamic matrix control - a computer control algorithm. In Proceedings of the Automatic Control Conference. San Francisco, CA, 1980 9. El Gahoui L, Scorletti G. Control of rational systems using linear-fractional representations and linear matrix inequalities. Automatica 1996; 32(9): 1273–1284 10. El Ghaoui L, Niculescu S (eds.) Advances in linear matrix inequality methods in control. Advances in Design and Control, SIAM, Philadelphia, PA, 2000 11. Fukuda M, Kojima M. Brach-and-Cut algorithms for the bilinear matrix inequality eigenvalue problem. Comput Optim Appl 2001; 19: 79–105 12. Gahinet P, Apkarian P. A linear matrix inequality approach to H1 control. Int J Robust Nonlinear Control 1994; 4: 421–448 13. Gahinet P, Apkarian P, Chilali M. Affine parameterdependent lyapunov funtions and real parametric uncertainty. IEEE Trans Autom Control 1996; 41(3): 436–442 14. Gahinet P, Nemirovski A, Laub AJ, Chilali M. LMI Control Toolbox. Technical report. The Mathworks Inc. 1995 15. Goodwin G, Graebe S, Salgado M. Control System Design. Prentice Hall, Upper Saddle River, NJ, 2001 16. Helmersson A. Parameter dependent lyapunov functions based on linear fractional transformation. In 14th IFAC World Congress. Beijing (China), 1999, pp. 537–542 17. Iwasaki T, Hara S. Well-posedness of feedback systems. IEEE Trans Autom Control 1998; 43(5): 619–630 18. Kothare KV, Balakrishnan V, Morari M. Robust constrained model predictive control using linear matrix inequalities. Automatica 1996; 32(10): 1361–1379 19. Kočvara M, Stingl M. A code for convex nonlinear and semidefinite programming. Optim Methods Soft 2003; 8(3): 317–333 20. Ljung L. System Identification. Theroy for the User, 2nd edn. Prentince-Hall, Englewood Cliffs, 1999 21. Martínez M, Senent J, Blasco FX. Generalized predictive control using genetic algorithms (GAGPC). Eng appl Artif Intell 1998; 11(3): 355–368 22. Megías D, Serrano J, Kuznetsov AG. A sistematic method to enhance the robustness of stabilising receding-horizon predictive controllers. In European Control Conference (ECC), Karlsruhe (Germany), 1999

463

GPC Robust Design with LMIs and/or BMIs

23. Packard A, Doyle J. The complex structured singular value. Automatica 1993; 29: 71–129 24. Rawlings JB. Tutorial overview of model predictive control. IEEE Control Syst Mag 2000; 20(4): 38–52 25. Salcedo JV. GPCs in state space for non linear systems. PhD thesis, Technical University of Valencia, 2005 26. Salcedo JV, Martínez M, Sanchis J, Blasco X. Design of GPC's in state space. Automatica 2002; 42(3–4): 159–167 27. Sanchez Peña RS, Sznaier M. Robust systems, theory and applications. John Wiley and Sons, New York, 1998 28. Scherer C, Gahinet P, Chilali M. Multiobjective outputfeedback control via LMI optimization. IEEE Trans Autom Control 1997; 42(7): 896–911 29. Scherer CW. LPV control and full block multipliers. Automatica 2001; 37: 361–375 30. Scherer CW. LMI relaxations in robust control. Eur J Control 2006; 12(1): 3–29 31. Scherer CW. Relaxations for robust linear matrix inequality problems with verifications for exactness. SIAM J Matrix Anal Appl 2006; 27(2): 365–395 32. Toker O, Özbay H. On the NP-hardness of solving bilinear matrix inequalities and simultaneous stabilization with static output feedback. In Proceedings of the American Control Conference, Seattle, WA (USA). 2525–2526, Vol. 4, 1995 33. Tuan H, Apkarian P. Low nonconvexity-rank bilinear matrix inequalities: Algorithms and applications in robust controller and structure designs. IEEE Trans Autom Control 2000; 45(11): 2111–2117 34. VanAntwerp JG, Braatz RD. A tutorial on linear and bilinear matrix inequalities. J Process Control 2000; 10: 363–385 35. Xi Y. New design method for discrete-time multi-variable predictive controllers. Int J Control 1989; 49(1): 45–56 36. Yoon T, Clarke D. Observer design in receding-horizon predictive control. Int J Control 1995; 61(1): 171–191 37. Zhou K, Doyle JC. Essentials of Robust Control. Prentice Hall, Englewood Cliffs, NJ, 1998

After some straightforward operations12 the following equivalent condition is obtained:

Appendix

if and only if, a symmetric multiplier V is present in such a way that:

A.1

 þ 0

ð40Þ R(x) > 0,

Starting from that robust stability condition (11) is obtained: T 

^ PðÞ 0 ^ 0 PðÞ

TA A

TB

0

0

TB

TC

8 0 TD  > < ~ IB  @ 0 > TA B : 0 0

1

~ GðÞ

T

C  Q 0

0 Q

 ~ <0 GðÞ

0 TD 0



I AðÞ

1 91 > = C~ TC B A > ; D 0

 ~ ^ , ^ 1 ~ 2 ~¼ diag , 8

ð42Þ The following step is to derive a finite number of LMIs using the full block S-procedure [10,29]: Lemma A.1. (Full block S-procedure) Suppose that P is a subspace of Rn , R 2 Rln is a full row rank matrix, and U Rkl is a compact set of matrices of full row rank. Then, the family of subspaces can be defined as follows: ~

P U ¼ P \ KerðURÞ ¼ fx 2 P : URx ¼ 0g ¼ fx 2 P : Rx 2 KerðUÞg ð43Þ is a fixed symindexed by U 2 U. Suppose N 2 R metric matrix and P 0 a fixed subspace of P with dimðP 0 Þ  k, Then: nn

8U 2 U : P U \ P 0 ¼ f0g and N < 0 on P 13 U ð44Þ

V > 0 on KerðUÞ

Sufficient conditions for robust stability with LFR uncertainty

I AðÞ



TA

8U 2 U : N þ RT VR < 0 on P and

is equivalent to matrix inequalities: Q(x)  S(x)R(x)  1S(x)T > 0.





 B C ~ ^ , ^  ¼ diag ,  @ TC A A ,

Schur complement

The matrix inequality:   QðxÞ SðxÞ >0 SðxÞT RðxÞ

A.2

~ ¼ GðÞ

ð45Þ

To apply this result to condition (42) the following selections have been made:  n

o ~ :  ~ 21 ~ IsizeðÞ  U¼ ~ ! 0rowsðÞcolsðQÞ IsizeðÞ 0 ~ ~ R¼ 0rowsðÞcolsðQÞ 0 IsizeðÞ ~ ~ 1 0 Q 0 0 C B 0 Q N¼@ A 0 02sizeðÞ ~

 <0

12

Operations between LFR expressions. See [9]. A matrix N < 0 is said to be negative (positive) definite over a subspace P U i.e. N < 0 ( > 0) on P U , when given a basis of P U collected in the columns of matrix BP U , then BTP U NBP U < 0ð> 0Þ. 13

ð41Þ

464

J.V. Salcedo and M.Martı´nez

0

TA B TA A B B B 0 B B 0 B P ¼ Span B B 0 B B T B C B @ TC A

TB 0

0 TB

IsizeðÞ ^ 0

0 IsizeðÞ ^

0 TD

0 0

0

TD

1 0 TA B C C C 0 C C 0 C C C IsizeðÞ C C 0 C C C TC B A

C 0 0 D |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} 0

T

TB

0

0

1

B 0 TB TA B C B C B C B IsizeðÞ C 0 0 ^ B C B 0 IsizeðÞ 0 C ^ B C P 0 ¼ Span B C B 0 C 0 I sizeðÞ B C B T 0 0 C D B C B C @ 0 TD TC B A 0 0 D |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} ð46Þ Applying full block S-procedure to (42) the following necessary and sufficient conditions of quadratic stability with LFR dependence are obtained: Q 0 0 Q 0 ~ 

! 0

V

ð47Þ

V

!T

IsizeðÞ ~

T<0 !

~ 

>0

IsizeðÞ ~

~ 21 ~ 8

 ~ ¼ Co X ~ 1 ; 

ð49Þ

In [10], the satisfaction of this LMI only in the finite set X ~ has been justified as a sufficient condition: ~ 

!T

IsizeðÞ ~

V

~  IsizeðÞ ~

ð53Þ

with fixed hermitian matrices Llj . This includes spectral norm bounds kl k 1, positive realness constraints l þ Tl  0 and mixtures thereof as the most important cases. Assuming this kind of uncertainty, condition (48) can be replaced by the sufficient conditions:  Diagðyl1 ,yl2 , . . . ,ylbl Þ 0 V¼

bl X

ylj Llj ,

>0

~ 2X~ 8 

ð50Þ

ð54Þ

yl1 , . . . ,ylbl 2

j¼1

which are convex. Diag(  ) is a block diagonal matrix. Repeated scalar blocks (real or complex) could be represented by: 

 j I : j 2 C, ð j

!

 ð j

1 ÞNmj

j 1

1 ÞN0j

j 1

 ¼ 0,



  0, m ¼ 1, . . . ,j

ð55Þ

with fixed hermitian matrices:

adding the LMI: V11 < 0

     l  0, j ¼ 1,...,bl l ¼ l 2 CFl Fl : Tl I Llj I

ð48Þ

The second condition, related to V multiplier, is based on an infinite number of LMIs. To acquire a finite ~ is given as the convex number it is assumed that set 1 hull of a finite number of matrices: n o ~ ~ ~ X ~ ¼  1 , . . . ,r

~ Remark A.1: In [31] it is observed that when the set 1 is a convex hull of a finite number of matrices, it includes general polytopes of real or complex full matrices, or of structured matrices with scalar elements that are confined to polytopes in real or complex spaces. ~ can be This idea reveals that the assumption on set 1 compatible with the uncertainty structure proposed by authors (6). Remark A.2: Other procedures to take into account cases when the general uncertainty structure (6) could not be represented as a convex hull, are described in [31]. For example, complex full blocks generally are described by:

T0

TT

where V11 is the (1,1)-block of symmetric  multiplier  V ~  partied conformably integrated with ~ : Isize ðÞ   V11 V12 ð52Þ V¼ VT12 V22

ð51Þ

 Nmj 2 C22 with det Nmj < 0 or Nmj ¼ 0

ð56Þ

465

GPC Robust Design with LMIs and/or BMIs

Specifically r þ 2 LMIs are present. If the a priori 1-norm bound,  is given, the unknowns are Q and V. However, these LMIs can be also used for obtaining the minimum value of .15  now being an additional unknown. This is a single optimization problem subject to (59).

Assuming this kind of uncertainty, condition (48) using hermitian matrices Ymj m ¼ 0, . . . ,j , can be replaced by the sufficient conditions:  DiagðY1j ,Y2j , . . . ,Yl j Þ 0 V¼

j X

ð57Þ

Nmj  Ymj

Remark A.4: The comments of Remarks A.1, A.2 and A.3 could be also applied for the 1-norm case.

m¼1

which are convex.  is the matrix tensor product. Remark A.3: In a general case of the uncertainty presented (6) it could be necessary to use the three kinds of ~ : convex hull, full blocks and blocks presented for 1 repeated scalar blocks. [31] shows how to simultaneously manage all of them.

A.4

Proof: Starting from from closed-loop transfer matrix (30): T0 ðzÞ ¼ Gcl ðzÞH1 ðzÞRðzÞ

^ (11) can be The positiveness condition of PðÞ removed under the assumption that nominal plant ( ¼ 0) is asymptotically stable14. This result is established in appendix A Proposition A.1. This condition can be confirmed by obtaining the eigenvalues of the uncertain state matrix, A(), (14) when  ¼ 0. A.3

3   A  C  BM B   BP 7 6 < ¼ 4 ðI  OÞM I  O ðI  OÞP 5, M I P 3 2      A  C  BM B   BU 7 6 H ¼ 4 ðI  OÞM I  O ðI  OÞU 5 M I U 2

Once again, taking the LFR dependence for the Lyapunov matrix (15) and LFR uncertainty for the plant:     A B! AðÞ B! ðÞ ¼ Cy Dy! Cy ðDÞ Dy! ðÞ ð58Þ   B 1 þ ðI  D Þ ð C D! Þ Dy

 ~  IsizeðÞ ~

!T V

~  IsizeðÞ ~

!

ð60Þ

Transfer matrices RðzÞ and HðzÞ display the following state space representation obtained from (28):

Sufficient conditions for 1-norm bounding

and applying full block S-procedure, a sufficient condition based on a finite number of LMIs, which guarantees that 1-norm is robustly less than , is obtained [25]: 0 1T  BC B C B C BC 0 1 B C Q 0 0 0 0 BC B C B C Q 0 0 0C B C B 0 BC B C B C B 0 0 0C 0  2 IcolsðB! Þ BC B C B C B C 0 0 IrowsðCy Þ 0 A B C @ 0 BC B C 0 0 0 0 V BC B C B C @A

Proof of Proposition 5.1

Unless pole/zero cancellations are present, T0 poles are closed-loop poles and H1 ðzÞRðzÞ poles. The poles of H1 ðzÞRðzÞ are obtained from state space representation. Firstly, H1 is obtained: 2 3 A  C 0 B 0 0 ðI  OÞ 5 H1 ¼ 4 1 1 ðUÞ1 ðUÞ M ðUÞ ð62Þ 0

TA

B TA A B B B 0 B B C B y B B 0 B B 0 B B B 0 B B T B C B @ TC A C

0

TB

0

TA B! I

0 0

TB 0

Dy! 0

0 IsizeðÞ ^

0 0

0

0

IsizeðÞ ^

0 0

0 TD

0 0

TB B! D!

0 0

TD 0

0

1

TA B C C C 0 C C Dy C C C 0 C C<0 0 C C C IsizeðÞ C C 0 C C C TC B A D

~ 2 X ~ , V11 < 0 > 0 8 

14 This result is also valid if the closed-loop is asymptotically stable for any fixed value of uncertainty

ð61Þ

Minimum 1-norm bound given the uncertainty set X ~ .

15

ð59Þ

466

J.V. Salcedo and M.Martı´nez

The H1 ðzÞRðzÞ product is:

2

A  C 6 0 6 0 H1 ðzÞRðzÞ ¼ 6 6 4 0 ðUÞ1 M

0 0 0 0 ðUÞ1

 BM ðI  OÞM  A  C  BM ðI  OÞM ðUÞ1 M

3  B BP ðI  OÞ ðI  OÞP 7 7  7 B   BP 7 ðI  OÞ ðI  OÞP 5 ðUÞ1 ðUÞ1 P

ð63Þ

This state space representation is not minimal, but this fact can be better observed making a linear state transformation:

2

A  C 6 0 6 0 H1 ðzÞRðzÞ ¼ 6 6 4 0 ðUÞ1 M

0 0 0 0 ðUÞ1

0 0  A  C  BM ðI  OÞM 0

The states corresponding to the last two row blocks are not observable, and, as such, these will be removed. Moreover, the states corresponding to the second row block are always zero, and consequently, these will also be removed. Finally, the minimal state space representation is: H1 ðzÞRðzÞ ¼



  C  A 1 ðUÞ M

 ðUÞ1 P

 ð65Þ

The state matrix coincides with observer state matrix, as such, the term H1 ðzÞRðzÞ takes the same poles as the observer. Consequently, T0(z) poles, unless pole/ zero cancellations are present, are closed-loop poles and observer poles. &

Other proofs ^ > 0) Let this be a discrete Proposition A.1: (PðÞ system with uncertainty whose state matrix is A(). ^ 0 exists for which the system is Suppose a fixed 0 ,  ^ 0 Þ > 0 exists, which verifies stable, that is to say, Pð the inequality: ^ 0 ÞAð0 Þ  Pð ^ 0Þ < 0 Að0 ÞT Pð

ð66Þ

0 0 B ðI  OÞ 0

3  7 0 7    BP 7 7 ðI  OÞP 5 ðUÞ1 P

ð64Þ

^ then Then, if this inequality is verified for every , , ^ ^ PðÞ > 0 for all . Proof: Initially, this result for continuous time systems will be used [13]: ^ > 0Þ Let this be a continuous Lemma A.2: ðPðÞ system with uncertainty whose state matrix is A(). ^ 0 for which the system is stable, Suppose a fixed 0 ,  ^ that is to say, Pð0 Þ > 0 exists, which verifies the inequality: ^ 0 Þ þ Pð ^ 0 ÞAð0 Þ < 0 Að0 ÞT Pð

ð67Þ

^ then Then, if this inequality is verified for every ,  ^ ^ PðÞ > 0 for all . The following step is to make a change of variable to transform discrete case to continuous one. Starting from the discrete time Lyapunov stability condition: ^ ^ <0  PðÞ AðÞT PðÞAðÞ

ð68Þ

^ < 0 occurs such that: then Zð,Þ ^ ^ ¼ Zð,Þ ^  PðÞ AðÞT PðÞAðÞ

ð69Þ

467

GPC Robust Design with LMIs and/or BMIs

This discrete Lyapunov equation can be transformed into the continuous case following a change in the variables: Ac ðÞ ¼ ðAðÞ þ IÞ1 ^ ¼ ðI  Ac ðÞÞZð,ÞðI  Ac ðÞÞT =2 ^ þ PðÞA ^ c ðÞ ¼ Zc ð,Þ ^ AT ðÞPðÞ c

ð70Þ

^ is negative verifying that Ac() is stable and Zc ð,Þ ^ ^ definite for all ,  if, and only if, A() and Zð,Þ satisfy the same conditions. From the Lyapunov continuous equation, it can be concluded that: ^ þ PðÞA ^ c ðÞ < 0 ATc ðÞPðÞ ^ > 0 for all . ^ this implies that PðÞ

^ 8,

ð71Þ &