LPV identification of a turbocharged diesel engine

Applied Numerical Mathematics 58 (2008) 1553–1571 www.elsevier.com/locate/apnum

LPV identification of a turbocharged diesel engine ✩ J.V. Salcedo ∗ , M. Martínez Department of Systems Engineering and Control, Universidad Politécnica de Valencia, Camino de Vera S/N, 46022 Valencia, Spain Available online 25 October 2007

Abstract This work proposes a methodology of identifying linear parameter varying (LPV) models for nonlinear systems. First, linear local models in some operating points, by applying standard identifications procedures for linear systems in time domain, are obtained. Next, a LPV model with linear fractional dependence (LFR) with respect to measured variables is fitted with the condition of containing all the linear models identified in previous step (differential inclusion). The fit is carried out using nonlinear least squares algorithms. Finally, this identification methodology will then be applied to a nonlinear turbocharged diesel engine. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Nonlinear systems; Identification; LPV models; Linear fractional transformation; Nonlinear curve fitting

1. Introduction The identification of local linear models of nonlinear systems can be carried out by means of the same identification techniques as those developed for linear systems [12]. Thus, in order to adequately represent a given nonlinear system in a specific area of operation it seems to be necessary to obtain by identification a sufficient number of linear models at different operating points in such a way that the possible gain variations, settling time, overshoot, etc., are reflected in the identified models. When all the models have been identified, the problem arises of how to integrate them in order to be able to correctly design either a single controller or several that guarantee certain specifications in the area under study. In this paper we propose to solve this problem by means of obtaining a linear model with time variable parameters (LPV) that integrates all the identified linear models at the different operating points. It must be noted that LPV models will have an LFR dependence with respect to time-varying parameters [7,18]. In the analyzed literature, the identification of LPV models have several approaches. Some works try to obtain directly a LPV model from simulated or experimental data. In [14,11,13,1,8] the mean square prediction error is minimized in order to fit the LPV model using the values of inputs, outputs and time-varying parameters, following the main ideas of prediction error method (PEM) used for linear time invariant systems [12]. In Refs. [14,13,1,8] the optimization problem is recast as a least square problem, whereas in [11] a nonlinear least square algorithm must be ✩

Partially supported by FEDER-CYCIT projects: DPI2004-08383-C03-02 and DPI2005-07835.

* Corresponding author.

E-mail address: [email protected] (J.V. Salcedo). URL: http://ctl-predictivo.upv.es (J.V. Salcedo). 0168-9274/$30.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2007.09.005

1554

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

Nomenclature diag(·) a block diagonal matrix whose entries are the matrices between parenthesis eq yj

the j th component of output vector at some operating point (equilibrium)

eq

ul

y˜ eq u˜ eq Co(·)

the j th component of input vector at some operating point the output vector at some operating point the input vector at some operating point the convex hull of a set

used. The proposed methodologies in [8,1,11] can be applied for any number of time-varying parameters, however the proposed ones in [14,13] are only valid of one time-varying parameter. On the other hand, in [20] a method derived from subspace identification methods [12] is performed in order to obtain a LPV in state space representation. The main limitation of this method is that it considers affine dependence of LPV model with respect to time-varying parameters. Also, [19] proposes an identification method for interpolating an LPV model composed by two terms: first term is a parametric model derived from a linear combination of a priori fixed LPV models, and the second term is a nonparametric model which, a priori, is ∞-norm bounded. The interpolation is carried out with LMIs using the Carathéodory–Féjer method. This last reference and [13] propose robust identification methods, in contrast with other works. Finally, in [21] a LPV model is identified in two steps. First, local linear frequency domain models are identified. And second, a LPV model in state space is fitted departing from the local models. This methodology is well developed only for a SISO system. The authors following the structure of last reference and the ideas of [2,3], present in this paper an alternative method for identifying LPV models, based in identified linear local models in state space domain for MIMO nonlinear systems. This methodology seems easier and systematic compared with the approach based on directly identify an LPV model [14,11,13,1,8,20,19], which requires, a priori, a previous effort to determine an adequate dependence with respect to time-varying parameters, and the fact that the optimization fitting must be applied to all the available data, which can be very expensive from a computational point of view. As opposite, in this paper, once the identified local models are available, a LPV model in state space with LFR dependence with respect to measured variables is fitted using nonlinear least squares, with the condition of fitting the identified local models. Besides, the values of systems signals at some operating points will be fitted using functions with LFR dependence. This idea is not performed in the literature analyzed. These values are quite important in order to ensure that the obtained LPV model will adequately reproduce signals values of the nonlinear system in the transients and in the steady states. Finally, this work presents an application for a turbocharged diesel engine, a highly nonlinear system. This paper is structured as follows: • • • • •

In nomenclature some definitions are presented. In Section 2 the LPV systems are introduced. In Section 3 the LPV identification method is presented. Section 4 develops the application for a nonlinear system: turbocharged diesel engine. And finally, Section 5 contains the conclusions.

2. LPV systems LPV systems are linear systems whose internal representation matrices are a function of one or various time-varying parameters λi (k): ˜ ˜ ˜ x(k + 1) = A(k)x(k) + B(k)u(k), y(k) = C(k)x(k) + D(k)u(k),         A(k) = f A Λ(k) , B(k) = f B Λ(k) , C(k) = f C Λ(k) , D(k) = f D Λ(k) ,   Λ(k) = diag λ1 (k)I s1 , λ2 (k)I s2 , . . . , λr (k)I sr .

(1)

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

1555

Isj is the identity matrix of size sj . The λi parameters can represent, for example, physical parameters or system output signals. For this work, a specific type of LPV system was used, being that in which the functional dependence of the f ∗ (f A , f B , f C , f D ) functions is linear-fractional (LFR) [6]:  −1 (2) f ∗ = M ∗ + L∗ Λ(k) I − Q∗ Λ(k) R ∗ this dependence allows the function so defined with respect to the parameters λi (k) to be rational. A joint LFR dependence is generally used for the four matrices:        −1 A B A(Λ(k)) B(Λ(k)) BΛ Λ(k) I − D Λ Λ(k) (C Λ D Λ,u ). = + (3) D y,Λ C(Λ(k)) D(Λ(k)) C D For example, if matrix A(Λ(k)) has order 1 and matrix Λ(k) is of order 2 with parameters λ1 (k) and λ2 (k) giving:     0 λ1 (k) a Λ(k) = a + (bΛ,1 bΛ,2 ) 0 λ2 (k)    −1  λ1 (k) cΛ,1 0 , × I 2 − dΛ cΛ,2 0 λ2 (k)   λ1 (k) λ2 (k) a Λ(k) = a + bΛ,1 · cΛ,1 + bΛ,2 · cΛ,2 . (4) 1 − dΛ λ1 (k) 1 − dΛ λ2 (k) The LFR dependence is reduced to the affine dependence if matrix D Λ is zero. In the foregoing example this is easily seen when dΛ is eliminated. Therefore, LFR dependence is more general than affine one. This LFR dependence is of general application since it can be used to give good representations of a large number of nonlinear systems [6,5,7], taking as λi (k) parameters specific system output signals involved in the nonlinear terms of their model. To ensure that the LPV model gives a good representation of the nonlinear system, differential inclusion techniques are employed [4]. These techniques are based on including the dynamic trajectories of the nonlinear system within those of the LPV model. This result shows that an LPV model with LFR dependence could adequately represent a nonlinear system for which any a priori model is available, although it is possible to carry out diverse identification experiments on it. The next step would then be to suggest a methodology to fit by identification this LPV model with LFR dependency. 3. Identification methodology of an LPV model 3.1. Literature review Refs. [2,3] present a methodology that consists of trying to fit an LPV model with LFR dependency that reproduces the different linearized models at each of the operating points. 3.1.1. Identification of local models For example, in [2] models are identified for a nuclear reactor at different operating points. These models are of the first order for the different outputs. Thus, for each output yj 1 the following linear models are obtained as a function ˜ of the input vector u(k) with the following structure:2 ˜ x(k + 1) = a(pi )x(k) + b(pi )u(k), ˜ + 1). yj (k + 1) = cj (pi )x(k + 1) + d j (pi )u(k

(5)

The parameter pi represents the value of the power of the reactor (p) at the operating point, which is a measurable signal. In all, there are nPF identified operating points. Given these models, a fitting is made to a linear model that presents LFR dependence with respect to the variable parameter λ(p) function of the power of the reactor: 1 The signals y , x and u ˜ represent the increase with respect to their operating point. j 2 The output equation is proposed at the instant k + 1 since the objective of the local model is to estimate the output at k + 1 knowing the operating

point at k, which depends on certain outputs of the system at k. Consequently, this local model (a, b, c and d) depends on instant k.

1556

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

   −1 aˆ λ(p) = a + bΛ λ(p) 1 − dΛ λ(p) cΛ ,  −1   bˆ λ(p) = b + bΛ λ(p) 1 − dΛ λ(p) d Λ,u ,  −1   cˆ λ(p) = c + dy,Λ λ(p) 1 − dΛ λ(p) cΛ ,  −1   dˆ λ(p) = d + dy,Λ λ(p) 1 − dΛ λ(p) d Λ,u .

(6)

This linear model with LFR dependence can be interpreted as a linear model variable with the operating point. This can be expressed more concisely (3):       ˆ    −1 a b a(λ(p)) ˆ b(λ(p)) bΛ λ(p) 1 − dΛ λ(p) (cΛ d Λ,u ). + G λ(p) = (7) = ˆ d c d y,Λ c(λ(p)) ˆ d(λ(p)) In the foregoing reference a, b, c, d, bΛ , dy,Λ , cΛ , d Λ,u , λ(p1 ), . . . , λ(pnPF ) are obtained by calculating the least squares solution of the following index: J=

nPF    2   2   2   2  a(pi ) − aˆ λ(pi ) + b(pi ) − bˆ λ(pi ) 2 + c(pi ) − cˆ λ(pi ) + d(pi ) − dˆ λ(pi ) 2 .

(8)

i=1

With the restriction that |λ(pi )| < 1 ∀pi . This problem of the least squares fitting is known in mathematical literature as NonLinear Least Squares Curve Fitting. The variable linear model with the operating point (6) is a particular case of a linear parameter varying system, since the power of the reactor is a system output signal, it therefore changes with time, p(k), when it passes from one operating point to another, causing parameter λ(p) to be also variable in time. Therefore, when this model (6) is used, what is really being used is the LPV model with LFR dependence:   λ p(k) = λ(k),  −1   aˆ λ(k) = a + bΛ λ(k) 1 − dΛ λ(k) cΛ ,  −1   bˆ λ(k) = b + bΛ λ(k) 1 − dΛ λ(k) d Λ,u ,  −1   cˆ λ(k) = c + dy,Λ λ(k) 1 − dΛ λ(k) cΛ ,  −1   (9) dˆ λ(k) = d + dy,Λ λ(k) 1 − dΛ λ(k) d Λ,u . The fitting of the linear model with LFR dependence to the local linear models could be improved, in theory, with a greater number of parameters to be fitted. The number of parameters depends on the complexity of the matrix Λ(k), comparing the proposed LPV model (9) with the general case (3) Λ(k) = λ(p(k)) is obtained. It would be possible to introduce a greater number of parameters in the fit by using a more complex matrix, such as, for example, Λ(k) = diag(λ1 (p(k))I s1 , . . . , λr (p(k))I sr ), with the objective of improving the initial fit if this one was not satisfactory. 3.1.2. Signal values at operating points It is important to emphasize that the value of λ(p) is only determined at the operating points analyzed. However, the interpolation of these values could provide the corresponding values to other operating points. Ref. [2] does not explain how the values of λ(p) are calculated or interpolated from λ(pi ) and λ(pi+1 ) when p is contained between the two operating points pi and pi+1 . Also, the works [2,3] does not make a fit for the systems signals in the operating points. Consequently, there is no guarantee that the fitted LPV model can really reach these operating points. This problem has been mentioned in the introduction. 3.2. Developed methodology 3.2.1. Operating points fitting We propose to fit by the least squares method these values at the operating points by means of matrix equations with LFR dependence with respect to λ(p), as has been done with the parameters of the linear model variable with operating point (6):

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

 −1 eq yˆj (pi ) = yj,0 + yj,1 λ(pi ) 1 − yj,2 λ(pi ) yj,3 ∀j, i,  −1 eq uˆ l (pi ) = ul,0 + ul,1 λ(pi ) 1 − ul,2 λ(pi ) ul,3 ∀l, i.

(10)

The parameters yj,0 , yj,1 , yj,2 , yj,3 , ul,0 , ul,1 , ul,2 and uj,3 are obtained optimizing this quadratic index:  nPF    eq 2  eq 2 eq eq yj (pi ) − yˆj (pi ) + uj (pi ) − uˆ j (pi ) J1 = i

j

eq yj (pi )

1557

(11)

l eq uj (pi )

where and are the real values of output and input signals at the operating points. With this additional fit an LPV model which matches the signals values at operating points will be obtained. This fit is not required by the methods that obtain directly an LPV model [14,11,13,1,8,20,19]. The operating point, q(k), for the general case will be defined by all the time varying parameters (1): T  (12) q(k) = λ1 (k), . . . , λr (k) . The different operating points analyzed will be represented by q i , which contains the values of time-varying parameters at each operating point. 3.2.2. Local models identification It is proposed to carry out the identification of local models at each operating point (q i ) for each of the outputs, yj , ˜ employing a model in internal representation. Subsequently, this separately, with respect to each of the inputs (u) model in internal representation will be transformed into the structure of the MISO model in state space used by [17]: ⎤ ⎡ 0 0 ... 0 −a0,j ⎢ 1 0 ... 0 −a1,j ⎥ ⎥ ⎢ ⎢ 0 1 ... 0 −a2,j ⎥ , Aj (q i ) = ⎢ ⎥ ⎥ ⎢ .. . . . . . .. ⎦ ⎣. . .. . . 0 ⎡ ⎢ ⎢ B j (q i ) = ⎢ ⎣ C j (q i ) = (0

...

...

1 −anj −1,j

nj ×nj

b0,j,1 b1,j,1 .. .

b0,j,2 b1,j,2 .. .

... ... .. .

b0,j,m b1,j,m .. .

bnj −1,j,1

bnj −1,j,2

...

bnj −1,j,m

···

0 1)1×nj .

⎤ ⎥ ⎥ ⎥ ⎦

,

nj ×m

(13)

Obviously C j matrices does not depend on operating point (q i ). The main advantage of this model is that it possesses a lower number of parameters than the identified internal representation model, which, in general, will possess some matrices Aj and C j full of coefficients. On grouping the models corresponding to all the outputs, the complete identified model is obtained at each operating point: ⎡ ⎤ ⎡ ⎤ A1 0 . . . 0 B1 ⎢ .. ⎥ ⎢ B2 ⎥ ⎢ 0 A2 . . . . ⎥ ⎢ ⎥ ⎥ , B(q ) = , A(q i ) = ⎢ ⎢ .. ⎥ i ⎢ . ⎥ . . ⎣ . . . . ⎦ ⎣ . . . 0 ⎦ B n r×m 0 . . . 0 An r×r ⎡ ⎤ C1 0 . . . 0 ⎢ .. ⎥ n  ⎢ 0 C2 . . . . ⎥ ⎢ ⎥ C=⎢ . , r = nj , (14) ⎥ .. .. ⎣ .. . . 0 ⎦ j =1 0 . . . 0 C n n×r ˜ x(k + 1) = A(q i )x(k) + B(q i )u(k), ˜ + 1) = Cx(k + 1). y(k

(15)

1558

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

This model could be nonminimal for some operating point. This case does not represent a problem since if there exist other operating points where the model is minimal, the LPV model to be fitted would have these number of states in order to accordingly represent all the local models. When the identified models have been obtained for all the operating points of interest, the least squares fitting is carried out to obtain the LPV model with LFR dependence with respect to the outputs; which is to say, the different measurable outputs are taken as λi (k):     ˜ , Λ(k) = diag y1 (k)I s1 , y2 (k)I s2 , . . . , yn (k)I sn  Λ y(k)     ˆ y(k) ˜ ˜ ˜ x(k + 1) = A x(k) + Bˆ y(k) u(k),   ˜ + 1) = Cˆ y(k) ˜ y(k x(k + 1).

(16)

(17)

It should be noted that the C matrices of all the linearized models are identical when the internal representation (15) ˆ y(k)) ˜ is utilized, which imply that in the output equation of the LPV model the matrix C( is constant and equal to C. To do the fit the following index will be optimized: J2 =

nPF     2   2  A(q i ) − A ˆ y(q ˜ i ) + B(q i ) − Bˆ y(q ˜ i) . 2

2

(18)

i=1

3.2.3. Global LPV model For a general nonlinear system, the fitted LPV model (17) and the fitted signals values at operating points (10) have the following expressions3 :     ˆ Λ(k) x(k) + Bˆ Λ(k) u(k), ˜ x(k + 1) = A ˜ + 1) = Cx(k + 1), y(k   −1 eq  yˆ˜ Λ(k) = y 0 + y 1 Λ(k) I − y 2 Λ(k) y 3 ,   −1 eq  uˆ˜ Λ(k) = u0 + u1 Λ(k) I − u2 Λ(k) u3 .

(19)

(20)

When the parameters of the LPV model and those of the operating points have been fitted, the nonlinear system can be represented using the ideas of [15] by the global LPV model:          eq  ˆ Λ(k) x(k) − x eq Λ(k) + Bˆ Λ(k) u(k) ˜ x(k + 1) = x eq Λ(k) + A − uˆ˜ Λ(k) ,     eq  ˜ + 1) = yˆ˜ Λ(k) + C x(k + 1) − x eq Λ(k) , y(k

(21)

where x eq (Λ(k)) is the value of the state of the LPV model in equilibrium, whose dependence with respect to Λ(k), in theory, is unknown [15]. From Eq. (21) it can be appreciated that for this global LPV model to adequately represent the identified nonlinear system, there would have to be some kind of relationship between the parameters of the LPV model and those of the operating points. In fact, in [15] it is indicated that the global LPV model correctly represents the nonlinear system if its linearized model at each operating point coincides with the corresponding identified local linear model. These authors indicate that this condition can be verified if x eq (Λ(k))4 satisfies a set of differential equations in partial derivatives. Once the equations in partial derivatives have been obtained, the existence of their solution will be analyzed, and when solutions exist a particular one will be determined. As will be seen later, the conditions that guarantee that the equations in partial derivatives have a solution establish that there must be coordination between the fit of the parameters of the LPV model and the fit of the operating points. 3 y˜ represents the output increment vector with respect to some operating point. 4 So far unknown.

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

1559

3.2.4. Global LPV model fitting The state of the model in equilibrium possesses, as we have already said, an unknown expression. Nevertheless, [15] propose obtaining it by imposing the condition that the linearized model of (21) around each operating point (y˜ eq , u˜ eq ) should coincide with the fitted linearized model (17). The linearized model of (21) is:   ˜ ∂x eq (y(k))  ˜ x(k + 1) − y(k)  ˜ ∂ y(k) ˜ y(k)= y˜ eq     ˜ ∂x eq (y(k)) eq  ˆ ˜ y(k) = A(y˜ ) x(k) −  ˜ ∂ y(k) ˜ y(k)= y˜ eq     ˜ ∂ u˜ eq (y(k))  ˆ y˜ eq ) u(k) ˜ ˜ + B( −  y(k) ,  ˜ ∂ y(k) ˜ y(k)= y˜ eq     ˜ ˜ ∂ y˜ eq (y(k)) ∂x eq (y(k))  ˜ + 1) − ˜ ˜ y(k = C x(k + 1) − y(k) . (22) |y(k)= ˜ y˜ eq y(k)  ˜ ˜ ∂ y(k) ∂ y(k) ˜ y(k)= y˜ eq    =1

For this model to coincide with (17) ∀y˜ eq it is necessary and sufficient that: eq eq  eq eq  ˆ y˜ eq ) ∂x (y˜ ) = B( ˆ y˜ eq ) ∂ u˜ (y˜ ) , I − A( ∂ y˜ ∂ y˜ ∂x eq (y˜ eq ) C =I ∂ y˜

or equally:    eq eq  ˆ y˜ eq )) ∂x eq (y˜ eq ) ˆ y˜ eq ) ∂ u˜ (y˜ ) (I − A( B( ∂ y˜ . = C I ∂ y˜

(23)

(24)

It is important to emphasize that this vectorial differential equation is in partial derivatives, since both the equilibrium state and the inputs in equilibrium are derived with respect to all the outputs. It is therefore possible to try to obtain the expression of the states in equilibrium for this vectorial equation in partial derivatives to be satisfied, since after eq ˆ y(k)), ˆ y(k)), ˜ ˜ ˜ the fitting of models and points of equilibrium A( B( C and uˆ˜ (y(k)) are obtained. ˆ y˜ eq ) and n rows of C, and r unknowns This vectorial equation is composed of r + n scalar equations, r rows of A( corresponding to each of the components of x eq . Therefore, there are more equations than unknowns. For the system to have a solution, it is necessary that at least n equations are linearly dependent on the other r for all y˜ eq . In Appendix A this is justified for the case of two outputs (n = 2), n equations are dependent on the rest if and only if:   1   ˆ ˜ eq )  ˆ eq eq n1 −1 rows B 1 (y ∂ u˜ (y˜ ) 1 1− i=0 aˆ i,1 (y˜ eq ) = ,  0 ∂y eq 1 ˆ ˜ ) rows B 2 (y   1   ˆ ˜ eq )  ˆ eq eq n2 −1 rows B 2 (y ∂ u˜ (y˜ ) 1 1− i=0 aˆ i,2 (y˜ eq ) = . (25)  0 ∂y eq 2 ˆ ˜ ( y ) B rows 1 These conditions are no more than two equations in partial derivatives where a relation is established between the parameters of the LPV model and those of the fitted operating points. If these sets of parameters are fitted independently, then it is almost certain that these equations will not be verified. To guarantee that they are complied with it will be necessary to impose them as a condition during the fitting phase, which implies that the fittings of the parameters of the LPV model and those of the operating points must be carried out simultaneously and in coordination imposing (25). These conditions, as can be appreciated, establish a relationship between the parameters of the LPV model and the derivatives of the of the operating points, i.e. they provide additional information to that which was available during the independent fitting, since it is now known what the values of the derivatives of the inputs at the operating points must be. Let us now return to the equations in partial derivatives that allow us to obtain x eq (24). In order to guarantee the existence of solutions, besides requiring that the number of equations coincides with the number of unknowns,

1560

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

since they are equations in partial derivatives, the Cauchy–Schwarz conditions must be satisfied in order to obtain a solution: ∂ 2 x eq (y˜ eq ) ∂ 2 x eq (y˜ eq ) = ∀k = j = 1, . . . , n, (26) ∂yj yk ∂yk yj which acquire the form:     eq eq  ˆ y˜ eq )) † B( ˆ y˜ eq ) ∂ u˜ (y˜ ) (I − A( ∂yk ∂ C Vk  ∂ =

∂yj ˆ y˜ eq )) (I − A( C

 V l = Vl (1)

···

† 

eq

ˆ y˜ eq ) ∂ u˜ (y˜ B( ∂yj Vj

∂yk T Vl (n) ;

eq

)

∀k = j = 1, . . . , n 

Vl (i) =



1, 0,

(27)

i = l, i = l,

where † indicates pseudoinverse matrix. Consequently, to guarantee the existence of x eq it is necessary that during the coordinated fitting phase of both sets of parameters conditions (25) and (27) are imposed. If we wanted to obtain the expression of x eq we would only have to integrate the equations in partial derivatives (24). This integration is solved in Appendix A. When we have ensured that the global LPV model adequately represents the identified local models, we confirm that all the dynamic trajectories of these models are included in the trajectories of the global LPV model. This will satisfy the condition of differential inclusion [4]. 3.2.5. Steps of the proposed method The proposed method has three main steps: (1) The signals values at the operating points are fitted using the cost index (11) with the LFR dependence (10). (2) The identified local models are fitted using the cost index (18) with the LPV model (17). (3) Finally, the global LPV model (21) is obtained using conditions (25) and (27). 3.2.6. Comparison to existing methods If we compare this new method with those ones which directly obtain an LPV model [14,11,13,1,8,20], we can highlight that: • In all methods (including the new one) an optimization algorithm is used to fit the model parameters, thus minimizing a certain cost index. This algorithm should be able to solve nonlinear least squares fitting problems. In the case of methods which directly obtain an LPV model, the mean square prediction error is minimized, whereas in the proposed method, two different indexes are minimized: one for fitting the signals values at the operating points, and another one for fitting the identified local models. • Those methods directly obtaining an LPV model require setting the structure of that given model beforehand. The method we propose here only requires setting the LFR dependence of matrix Λ(k) with respect to time-varying parameters. • The fact that all the available experimental data are used for fitting in the methods which directly obtain an LPV model implies a high computational cost for the optimization algorithm. However, in the proposed method we only need to use for fitting, on the one hand, the signals values at the operating points and, on the other hand, the parameters of the identified local models. • Moreover, the proposed methodology has the advantage of adequately representing all the areas of the nonlinear system operating region only by including the appropriate number of operating points. It is more complex to achieve this with other methods since, a priori, it is not easy to determine how many experimental input/output data we need to obtain in each area so as to get a good representation.

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

1561

3.3. Choice of Λ(k) matrix structure The Λ(k) (16) matrix normally depends on all the measurable output signals. However, no indication has been given of how to choose or decide the sj sizes of the identity matrices associated with each of them. Clearly, the greater the size of these matrices, the greater the complexity of the LFR dependence. We propose a practical method of choosing the sj sizes. Initially, they are all taken to be equal to 1, and we proceed with the solution of the optimization problem associated with the fitting. If the fit thus obtained is satisfactory the process ends here. If not, some of the sizes are increased and the fitting process is repeated. This procedure can be repeated until a satisfactory result is obtained. In this work a satisfactory result is obtained if the minimum value of J2 (18), minimum square error, provides a mean square error by number of fitted parameters, ep , lesser or equal to a tolerance  > 0 a priori fixed: √ J ep  < (28) npar being npar the number of parameters to be fitted. Two fittings are performed in the process of obtaining the global LPV model: one for the LPV model parameters and another for those of the operating points. When the above proposed method is applied in both fittings, using the same structure for the Λ(k) matrix is not a necessary requirement for a satisfactory result. Therefore, from this point on, ΛM (k) will represent Λ(k) matrix for the fitted LPV model parameters and Λeq (k) the corresponding to the fitted parameters pertaining to the operating points. Moreover, for the operating points parameters, in order to analyze if the fit is satisfactory, the a priori tolerance pf is fixed. 3.4. Validation of the global LPV model If the global LPV model obtained with the previous methodology will be used for designing controllers for the nonlinear system, in order to guarantee an adequate controller performance, it will necessary to develop a method to validate the global LPV model. The validation method must guarantee that the signal values at the operating points used in the fitting are accurately reproduced by the operating point fitting, and the linear local models around these operating point are quite similar to the ones provided by the fitted LPV model. Also, in other operating points different from used in fitting, both signal values and linear local models must be adequately approximated by the global LPV model. Nevertheless, it must be expected that higher distance from operating points used in fitting worse reproduction will be achieved. If in that operating points the interpolated signals values and/or linear local models are poorly reproduced, that operating points will must be included in the fitting. The comparison between real signal values and the fitted ones is carried out computing their relative errors. If these errors are small a good fit is achieved. This method will be applied later in the example. The comparison between identified local models and local fitted models is carried out computing the relative error regarding the steady state gains, and representing in a same graph the frequency response of both models. Both tasks will be applied later in the example. 4. Application to a turbocharged diesel engine 4.1. Introduction Nowadays, turbocharged diesel engines produce more power and have better performance than a non-turbocharged diesel engine. However, due to the extra air consumed and higher temperatures and pressures inside the engine, they have the disadvantage of producing oxides of nitrogen (NOx ) in the exhaust system. The present regulations on the emission of these gases are quite strict, and manufacturers have opted to include in the design the recirculation of inert gases (H2 O and C2 O) between the exhaust manifold and the inlet manifold, the so-called EGR circuit. The recirculated air mass flow rate of these gases is designated as m ˙ EGR . It can be modified by means of the EGR valve in the EGR circuit. The main problem associated with this device is that it significantly reduces engine performance.

1562

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

Fig. 1. Turbocharged diesel engine diagram.

Additionally, in the near future the regulations will become much more strict, with the probable result that the present recirculation system will not be able to meet the new requirements regarding emission levels of oxides of nitrogen. In order to counteract the loss of power or torque due to the recirculation of gases, the turbine valve position can be modified (variable geometry turbines, known as VGT) so as to increase or reduce the pressure of the gases in the exhaust manifold. As a result, the power produced by the exhaust system can be varied so as to provide more or less power for the compressor. As the compressor power is increased the air mass flow also increases m ˙ a (which allows more fuel to be consumed) and as a consequence the engine produces more torque. The operating point defined by the manufacturer for the turbocharged diesel engine is controlled by the following variables: • Quantity of fuel injected (m ˙ f ), controlled by the accelerator pedal. • Engine speed (N ), controlled by engine torque and road resistant torque. For each operating point the manufacturer has a map that shows the airflow quantities (m ˙ a ) and inlet manifold pressures (Pa ) required to avoid emissions of smoke and oxides of nitrogen and also to guarantee a specific engine torque. Both values must be used as the references to determine the EGR and VGT valve positions for whatever control algorithm is implemented, thus achieving a reduction in the emission of contaminating gases and at the same time giving the level of power required by the driver. 4.2. Model of engine In cooperation with the Department of Internal Combustion Engines of the Polytechnic University of Valencia we have developed the first stage of what we hope will be a simplified model of the compressor-engine-circuit group of EGR+turbine (see Fig. 1), following a similar method as presented in [9,10]. The model contains both differential equations and static equations as well as the tables obtained from experiments carried out on the actual engine. Given its particular structure, we chose to implement it in the Simulink environment of MatLab with the objective of obtaining a model easy to simulate [16]. In performing the simulations certain variables can be modified according to requirements: (1) (2) (3) (4)

Engine velocity (N ). Fuel flow (m ˙ f ). % EGR valve position. % VGT valve position.

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

1563

Fig. 2. Operating zone of turbocharged diesel engine.

Of all the output variables provided by the Simulink model, the following are worthy of mention: (1) Air mass flow rate through compressor (m ˙ a ). (2) Inlet manifold pressure (Pa ). These are the two controlled outputs of the system. The model was approved by the staff of the Internal Combustion Engine Department of the University, who confirmed that the model performed in a very similar way to a real engine. 4.3. Identification of an LPV model Even if a Simulink engine model is available, it cannot be used to design controllers that attempt to stabilize the system in a given operating zone. For this we propose to identify an LPV model for turbo diesel engine from the Simulink model, by means of the methodology described in Section 3. The main difficulty with this methodology lies in the least squares fitting, since the structure of the Λ(k) matrix that the LPV model will have is generally unknown and the optimization problem will almost certainly be non-convex. To obtain the LPV model, the following engine speed and fuel flow values were selected: N = 1500 rpm and m ˙ f = 0.003 kg/s, corresponding to actual engine values. With this we now have a 2 × 2 system: control of VGT and EGR positions, controlled outputs m ˙ a and Pa . Using the Simulink model, we establish that when selecting N and mf , ma can vary between 90 and 120 kg/hour and Pa can vary between 1 and 1.2 bar. The operating points analyzed are shown in Fig. 2. It is important to point out that within the attainable ranges of each controlled output, it was not possible to obtain all their possible combinations. Nevertheless, the points shown in Fig. 2 are representative of the combinations obtainable. A higher number of attainable points could be used, but after the experiments it seemed reasonable to use only those shown. It must be remembered that the more points used, the greater the number of models, and therefore more complicated will be the process of fitting the LPV model. 4.3.1. Operating points fitting The values of the control variables (%EGR and %VGT) were obtained at the operating points shown in Fig. 2. Two PIDs were used in the calculation to stabilize the engine at the desired operating points. The PIDs were designed from a linearized model of the engine for the operating point Pa = 1.2 bar and m ˙ a = 100 kg/h. The first PID stabilizes Pa through %VGT, and the second m ˙ a through %EGR. This procedure is one of those habitually used to control multivariable processes. The values of the controlled variables, %VGT and %EGR, for the different equilibrium points were fitted by LFR dependence with respect to Pa and m ˙ a . To perform this fit the Λeq (k) matrix structure must be selected, in line with

1564

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

the methodology proposed in Section 3.3. After a series of attempts to fit %EGR and %VGT we observed that it was necessary to use the following structure for Λeq (k) (pf = 0.01) in order to ensure a good fit:   0 (Pa (k) − 1.1084)I 3 . Λeq (k) = 0 (m ˙ a (k) − 102.0147)I 2 The values 1.1084 and 102.0147 represent the mean values of Pa and m ˙ a , respectively, at the analyzed operating points. Their use made the optimization problem easier. 4.3.2. Local models identification and fitting The identification of the linearized models at the different operating points was performed in the following stages: (1) Generation of exciting input signals for EGR and VGT positions. It was decided to use pseudo random binary signals (PRBS). (2) Simulation of the Simulink model subjected to PRBS signals. (3) Identification of second order models in internal representation for both outputs. This selection gave a good fit. These identified models are transformed into the internal representation given in (13). As has already been stated, on using this particular representation, all the C matrices are identical and consequently it is not necessary to fit the ˜ matrix C(y(k)) since it coincides with C. Given the models in this internal representation, the final stage consists of fitting the local A and B matrices to an LPV model with LFR dependence by means of least squares optimization. This task was performed with the help of the MatLab Optimization Toolbox, which implements the nonlinear least squares fitting algorithm proposed by Levenburg and Marquardt. To ensure a good fit ( = 0.01) the ΛM (k) matrix structure must be selected, following the method described in Section 3.3. After several attempts the following fit was achieved:   0 (Pa (k) − 1.1084)I 4 . (29) ΛM (k) = 0 (m ˙ a (k) − 102.0147)I 4 The fit was performed, firstly, to obtain an LPV model for Pa , and, secondly, to obtain an LPV model for m ˙ a . The ˙ a in Eq. (B.1). Both equations are LPV model corresponding to Pa is shown in Eq. (B.1), and that corresponding to m placed at Appendix B. 4.3.3. Global LPV model fitting In Section 3 it was shown that the parameters of the LPV model and those of the operating points must satisfy certain conditions (25), (27). Thus, a re-fit is made of the operating point parameters, previously fitted individually, so that they comply with the first condition (25). This re-fit consists of adding to the fit index (11) additional terms that penalize quadratically the noncompliance with (25):       1 1     2 2 eq eq n2 −1 rows B 1 ∂ u rows B 2 ∂ u ˜ ˜ 1 −1 a 1 1 1−ni=0 1− a i,1 i,2 i=0 − − (30) + .   0 0 ∂y1 ∂y2 B2 B1 2

rows

2

rows

After this re-fit it became possible to verify approximately this condition, obtaining the following results:    −1 %EGR Λeq (k) = %EGR0 + %EGR 1 Λeq (k) I − %EGR 2 Λeq (k) %EGR 3 = 58.8642 + (47.9760 − 1.2211 6.0919 − 2.5747 − 4.8634) · Λeq (k) ⎛ ⎛ ⎞−1 ⎛ ⎞ ⎞ −3.9337 −19.2089 −13.8247 0.0222 −1.7026 17.3748 10.7303 −8.7396 −0.4030 −1.4471 ⎟ ⎜ ⎜ 9.2929 ⎟ ⎜ 5.3307 ⎟ ⎜ ⎜ ⎟ ⎜ ⎟ ⎟ × ⎜I 5 − ⎜ 0.1170 25.0404 4.1986 0.0874 0.8107 ⎟ Λeq (k)⎟ · ⎜ −0.9943 ⎟ , ⎝ ⎝ ⎠ ⎝ ⎠ ⎠ −1.4355 −2.2320 −1.2865 0.0303 −0.2015 3.5328 0.6392 0.8670 0.4841 −0.0138 0.0564 0.1271    −1 %VGT Λeq (k) = %VGT 0 + %V GT 1 Λeq (k) I − %V GT 2 Λeq (k) %V GT 3 = 50.6677 + (19.9639

− 5.4646

− 6.5723 2.0860

− 2.6115) · Λeq (k)

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

1565

Table 1 Relative errors between real and fitted values of %EGR Type

Pa

m ˙a

Rel. e. (%)

min max mean

1.05 1.2

97 115

0.0023 0.6268 0.0955

Table 2 Real and fitted values of %VGT and their relative error

⎛

Type

Pa

m ˙a

Rel. e. (%)

min max mean

1.1 1.15

99 99

0.0003 0.2765 0.0695

⎛

−8.0401 ⎜ ⎜ −6.0381 ⎜ ⎜ × ⎜I 5 − ⎜ −1.8889 ⎝ ⎝ −0.1459 −0.7276

⎞ ⎞−1 ⎛ ⎞ −7.2450 2.9884 −2.4865 −1.0396 −18.3501 −9.7158 −10.9696 −0.1663 0.6733 ⎟ ⎟ ⎜ 7.9306 ⎟ ⎟ ⎟ ⎜ ⎟ 7.0807 −8.1916 0.6268 0.4155 ⎟ Λeq (k)⎟ · ⎜ −1.2427 ⎟ . ⎠ ⎠ ⎝ ⎠ −0.2987 −1.8186 −0.1126 −0.0091 0.5182 −2.5912 1.2171 −0.1045 0.0103 0.2686

The other condition (27) was also approximately satisfied after the re-fit. 4.4. Validation of LPV model In this section a validation of previous LPV identified model is made, following the ideas presented in Section 3.4. The procedure is decomposed in 4 steps, 2 devoted to checking the signals fitting, and the other 2 related to analyze local models fitting. 4.4.1. Fitted values for input signals %EGR and %VGT at the operating points used in the fitting In Table 1 the maximum, the minimum and the mean value of the relative error between real and fitted values of input %EGR at some operating points used for fitting are showed. In Table 2 the same is made but for input %TGV. In all the cases the relative error is lesser that 1%. 4.4.2. Local models fitted at the operating points used in the fitting In Table 3 the maximum, the minimum and the mean the relative error (%) between static gains of the identified local models for outputs Pa and m ˙ a and the local models provided by the fitted LPV model5 are presented. At operating point (1.06, 95) the biggest relative errors are found: ≈ 1% for models corresponding to output Pa and ≈ 5% for models corresponding to output m ˙ a. In order to compare the dynamic behavior between both sets of models6 , the frequency response of larger singular value of transfer matrices for outputs Pa and m ˙ a corresponding to operating point (1.06, 95) are represented, which has the worst fit for the static gains. In Figs. 3 and 4 the frequency response of larger singular value corresponding to the identified local model and to the local fitted model are represented, for outputs Pa and m ˙ a respectively. As figures 5 That relative error is computed as the 2-norm of the difference between both gains divided by the 2-norm of static gain of identified local model and multiplied by 100. 6 Identified and fitted.

1566

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

Table 3 Relative error in 2-norm between static gains of identified local models and local ˙a fitted models for outputs Pa and m Type

Pa

m ˙a

Rel. e. gain Pa (%)

max min mean

1.05 1.2

93 100

1.0221 0.0296 0.3017

Type

Pa

m ˙a

Rel. e. gain m ˙ a (%)

max min mean

1.06 1.1

95 104.5

5.2629 0.0539 2.1136

Fig. 3. Frequency response of larger singular value for identified local model (continuous line) and for local fitted model (dashed line) for output Pa at operating point (1.06,95).

Fig. 4. Frequency response of larger singular value for identified local model (continuous line) and for local fitted model (dashed line) for output m ˙ a at operating point (1.06,95).

Table 4 Relative errors between real and fitted values of input signals %EGR and %VGT Type

Pa

m ˙a

Rel. e. %EGR (%)

min max mean

1.08 1.19

101 99

0.0223 4.6702 1.8504

Type

Pa

m ˙a

Rel. e. %VGT (%)

min max mean

1.19 1.11

99 105

0.0239 0.9851 0.4214

show, both responses are quite similar which guarantees the approximation between identified local models and local fitted models. Similar figures are obtained for the other operating points used in the fitting. 4.4.3. Fitted values of input signals %EGR and %VGT at operating points not used in the fitting In Table 4, for some operating points not used in the fitting, the maximum, the minimum and the mean value of the relative error between real and fitted values for input signals %EGR and %VGT are showed. It is appreciated that the relative errors are small, with exception of some operating points in which the errors reach values nearly to 5%. 4.4.4. Local models fitted at the operating points not used in the fitting The operating point that have the worst fit in the previous section is (1.19, 99), since it has a relative error of 5%, approximately, in the input signal %EGR. For such point in Figs. 5 and 6 frequency responses of larger singular values ˙ a respectively. corresponding to identified local models and fitted local models are represented, for outputs Pa and m

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

Fig. 5. Frequency response of larger singular value for identified local model (continuous line) and for local fitted model (dashed line) for output Pa at operating point (1.19, 99).

1567

Fig. 6. Frequency response of larger singular value for identified local model (continuous line) and for local fitted model (dashed line) for output m ˙ a at operating point (1.19, 99).

As figures show, the bigger discrepancy corresponds to m ˙ a model, in particular due to different values for medium frequencies, whereas at low/high frequencies both models are more similar. For other operating points the results are also similar. Concluding, along the previous lines it has been justified that the global LPV obtained for the turbocharged diesel engine provides a reasonable approximation of nonlinear system, better when the system stays near an operating point used in the fitting. The distancing from such point supposes a worse approximation, as it would be expected. 5. Conclusions • In the case of nonlinear systems a methodology was proposed to identify LPV models with fractional linear dependence whose dynamic trajectories contain those of the nonlinear process. • This methodology consists of obtaining diverse identified linear models at different operating points and subsequently by means of least squares fitting using the Levenberg and Marquardt algorithm obtaining an LPV model that attempts to reproduce the identified models at each operating point. • The aspect of the integration of the identified linear models with the values of the signals at each operating point was dealt with. A methodology was proposed based on two systems of equations in partial derivatives that ensure the agreement of both during the least squares fitting phase. • As an example, the proposed methodology was applied to a turbocharged diesel engine to demonstrate its practical utility. • The identified LPV model can be used for the design of nonlinear controllers or linear controllers with timevarying parameters. Acknowledgements We would like to thank the R&D+i Linguistic Assistance Office at the Universidad Politécnica de Valencia for their help in translating this paper. Appendix A. Proofs The vectorial differential equation (24) is a system of n + r equations with r unknowns. This system will have a solution if there are at least n equations dependent on the rest, or similarly:    eq eq  (I − A(y˜ eq )) (I − A(y˜ eq )) B(y˜ eq ) d u˜ d(y˜y˜ ) rank = rank . C C I

1568

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

˜ and C, the question lies in justifying For the case of two outputs, given the particular structure of the matrices A(y) that in the following matrix there are n = 2 rows dependent on the rest: eq eq   (I − A(y˜ eq )) B(y˜ eq ) d u˜ d(y˜y˜ ) ⎡ ⎢ ⎢ =⎢ ⎣

C I n1 − A1 (y˜ eq )

I 0

eq

u˜ B 1 (y˜ eq ) ∂∂y 1

eq

u˜ B 1 (y˜ eq ) ∂∂y 2

⎤

eq ∂ u˜ eq ⎥ u˜ eq ˜ 0 I n2 − A2 (y˜ eq ) B 2 (y˜ eq ) ∂∂y B ( y ) ∂y2 ⎥ 2 ⎥. 1 ⎦ 0...0 1 0 1 0 0 0...0 1 0 1 For the case in which the identified models have no integrators, it is then easy to justify that there are exactly two rows dependent on the rest if and only if:    1 ˜ eq ) ∂ u˜ eq (y˜ eq )  1  n1 −1 rows B 1 (y 1− i=0 ai,1 (y˜ eq ) = ,  0 eq ∂y1 ˜ B ( y ) 2 rows    1 ˜ eq ) ∂ u˜ eq (y˜ eq )  1  n2 −1 rows B 2 (y 1− i=0 ai,2 (y˜ eq ) = .  0 eq ∂y2 ˜ B ( y ) 1 rows These conditions can easily be extended to the case of the greater number of outputs. For the particular case in which one of the models corresponding to an output had an integrator, e.g. the number 1, then its corresponding condition would be:     B 1 (y˜ eq ) ∂ u˜ eq (y˜ eq ) 0 rows = . eq 0 ∂y ˜ B ( y ) 1 2 rows

To ensure that the system of equations in partial derivatives has a solution, the Cauchy–Schwarz conditions relative to the second derivatives must also be complied with. Once it is justified that the solution exists, the equations must be integrated. The following is an outline of how to obtain a solution in the case of having only two outputs, n = 2, although the method is perfectly viable for a greater number. Let the system of equations in partial derivatives be: eq ⎛ ˜ ⎞  † B(y) ˜ ∂ u˜∂y1(y) ∂x eq ˜ (I − A(y)) ⎝ ⎠ = f 1 (y1 , y2 ), = 1 C ∂y1 0 eq ⎛ ˜ ⎞   ˜ ∂ u˜∂y2(y) † B(y) ∂x eq ˜ (I − A(y)) ⎝ ⎠ = f 2 (y1 , y2 ). = 0 C ∂y2 1 It is important to observe that the functions f 1 and f 2 present an LFR dependence with respect to y1 and y2 . If the first equation is integrated with respect to y1 we obtain as a result that the state of equilibrium must have the form: # x eq = f 1 (y1 , y2 ) dy1 + g(y2 ) where g can be any arbitrary function of y2 . If this formula is substituted in the second differential equation: # dg ∂ = f2 (y1 , y2 ) f 1 (y1 , y2 ) dy1 + ∂y2 dy2 which permits function g to be obtained integrating with respect to y2 this last equation:  # #  # ∂ f 1 (y1 , y2 ) dy1 dy2 . g(y2 ) = f 2 (y1 , y2 ) dy2 − ∂y2 With which a solution of the system of differential equations in partial derivatives is:  # # #  # ∂ eq x = f 1 (y1 , y2 ) dy1 + f 2 (y1 , y2 ) dy2 − f 1 (y1 , y2 ) dy1 dy2 . ∂y2

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

1569

It is possible to justify that all the solutions to this system of differential equations in partial derivatives are differentiated by an arbitrary constant. The formula presented provides one of them and any other can be obtained by adding any constant. From this result it can be deduced that to solve this system of differential equations in partial derivatives it is only necessary to be able to integrate and derive matrix functions with LFR dependence. Derivation is the simplest of procedures, and the following formula can be deduced for this purpose:    −1 G Λ(z) = M + LΛ(z) I − DΛ(z) R,    dΛ(z) dG dΛ(z) Λ(z) 0 =L R+ L+L D 0 Λ(z) dz dz dz −1  dΛ(z)    dΛ(z)   D D dz D D dz R Λ(z) 0 × I− 0 Λ(z) 0 D R as G has a linear-fractional dependence with respect to z it is then assumed that Λ(z) only depends affinelly on z, which guarantees that the derivative once again has an LFR dependence with respect to z. The integration is slightly more complicated and its principal characteristic is that, in general, it does not lead to a representation with LFR dependence: # #     −1  G Λ(z) dz = M + LΛ(z) I − DΛ(z) R dz. If the D matrix is invertible it can then be shown that the integral has the expression: # #      −1   G Λ(z) dz = M + L −D −1 + D −1 I − DΛ(z) R dz # −1   −1  = (M − LD −1 R)z + D I − DΛ(z) R dz   = (M − LD −1 R)z − D −2 logm I − DΛ(z) R. Where the last step is based on the primitive of −D(I − DΛ(z))−1 being the matrix logarithm function logm(I − DΛ(z)). In the operations it is assumed that Λ(z) = I z. In the case of D not being invertible, i.e. having some zero eigenvalues, then the question becomes more complicated. In the first place, the real Schur decomposition must be calculated, which can be obtained with the schur function of MatLab, and order the elements of the matrix7 in such a way that the following structure is adopted:   D I D 1,2 T T D=T 0 D0     T  D I D 1,2 T 1f T . = T 1c T 1c 0 D0 T 1f D 0 being a matrix that only has zero eigenvalues, i.e. a nilpotent matrix, and D I an invertible matrix. With this decomposition, the following is obtained: −1   −1 I − D I z −D 1,2 z =TTz T Λ(z) I − DΛ(z) 0 I − D0z = T T1c z(I − D I z)−1 T 1f + T T1c z(I − D I z)−1 D 1,2 z(I − D 0 z)−1 T 2f + T T2c z(I − D 0 z)−1 T 2f . As the matrix D 0 is nilpotent then k > 0 exists, so that D k0 = 0 and Dk+1 = 0. In such a case: 0 (I − D 0 z)−1 = I + D 0 z + D 20 z2 + · · · + D k0 zk . After all these calculations we obtain: 7 For example, with the ordschur function available in MatLab 7.0.1 and later versions.

1570

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

G(z) = M + LT T1c z(I − D I z)−1 T 1f R + LT T1c (I − D I z)−1 D 1,2 (I z2 + D 0 z3 + · · · + D k0 zk+2 )T 2f R + LT T2c (I z + D 0 z2 + D 20 z3 + · · · + D k0 zk+1 )T 2f R as D I is now invertible, the method described for the case in which D is invertible can be applied. Additionally, it must take into consideration the following relationship: −j

−j

j −1 j −2 − D −2 − · · · − D I + D I (I − D I z)−1 . zj (I − D I z)−1 = −D −1 I z I z

Appendix B. Equations 

    APa ΛM (k) B Pa ΛM (k)

 −1 (C Pa ,Λ |DPa ,Λ,u ) = (APa |B Pa ) + B Pa ,Λ ΛM (k) I − D Pa ,Λ ΛM (k)   −6 1.4307 × 10−5 0 −0.8932 −8.3889 × 10 = 1.0000 1.8930 8.3613 × 10−6 −1.5091 × 10−5   3.4535 3.7696 −1.4358 −0.0313 0.1234 −0.0310 0.3511 −0.0053 + ΛM (k) −1.0673 −8.1126 −0.5341 0.1714 −0.1442 0.0234 −0.3481 0.0077 ⎛ 0.9013 ⎞−1 ⎛ 1.1135 −0.0316 −3.7984 0.4369 0.3318 −0.0657 0.0601 ⎞ 0.0667 0.0689 ⎟ ⎜ −0.1327 −0.4879 −0.5682 −1.9180 −0.1387 0.0161 ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ 3.0250 −0.9790 1.0172 0.7536 0.3338 −0.2098 −0.0734 ⎟ ⎜ 2.5561 ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ ⎜ −1.8874 1.1263 −27.7903 2.2520 −1.1746 −0.2564 −0.1860 0.1148 ⎟ ⎟ ⎜ × ⎜I 8 − ⎜ ⎟ ΛM (k)⎟ 0.2080 −1.5841 0.8877 0.1787 0.0204 −0.0295 −0.0039 ⎟ ⎜ 0.1788 ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ 2.8344 −1.4796 −0.5755 −0.1546 0.1036 0.0357 ⎟ ⎜ 1.0354 −2.8637 ⎟ ⎜ ⎝ ⎠ ⎠ ⎝ 0.6529 −0.9767 2.9480 0.6192 −0.3901 −0.1619 0.0703 0.0223 −2.1004 −3.4444 0.2053 −1.4653 −0.3456 −0.3443 0.1072 0.1085 ⎛ 0 0.0813 −3.901 × 10−6 3.9212 × 10−5 ⎞ ⎜ 0 0.0509 −3.8478 × 10−7 4.4156 × 10−5 ⎟ ⎜ ⎟ ⎜ 0 −0.0181 −3.4333 × 10−6 −4.7734 × 10−5 ⎟ ⎜ ⎟ ⎜ 0 −0.0508 6.3389 × 10−6 −1.8166 × 10−5 ⎟ ×⎜ ⎟, −6 −6 −9.8461 × 10 ⎟ ⎜ 0 0.0069 −3.7612 × 10 ⎜ ⎟ −6 −5 2.1101 × 10 ⎜ 0 −0.0176 6.1653 × 10 ⎟ ⎝ ⎠ −6 −6 0 −0.0129 1.6895 × 10 −2.2027 × 10 0 0.0285 −1.3909 × 10−5 2.1013 × 10−5     −1   Am˙ a ΛM (k) B m˙ a ΛM (k) = (Am˙ a |B m˙ a ) + B m˙ a ,Λ ΛM (k) I − D m˙ a ,Λ ΛM (k) (C m˙ a ,Λ |Dm˙ a ,Λ,u )   0 −0.87808 0.010696 −0.0095768 = 1 1.8778 −0.010729 0.0094919   1.4193 −1.9467 3.1852 −0.0114 −0.1242 0.1740 0.1909 −0.0092 + ΛM (k) −1.4103 1.9801 −3.1967 −0.0158 0.1235 −0.1731 −0.1909 0.0095 ⎛ 7.3291 ⎞−1 ⎛ 68.7948 −12.0761 0.4713 0.2603 0.4315 0.9797 −0.9592 ⎞ ⎜ −7.1668 −2.3819 −0.8765 ⎜ ⎜ ⎜ 12.3613 3.6694 ⎜ 2.9832 ⎜ ⎜ ⎜ 0.0577 ⎜ −0.0618 0.2964 ⎜ × ⎜I 8 − ⎜ 0.5241 ⎜ −0.5305 −0.9890 ⎜ ⎜ ⎜ 0.3927 ⎜ −0.3843 −1.1497 ⎜ ⎝ ⎝ −1.3418 −0.5501 0.1769 −0.0141 0.3142 −0.3245 ⎛0 0.242 0.0070 −0.1843 ⎞ ⎜ 0 −0.1706 −0.0058 0.0234 ⎟ ⎜ ⎟ ⎜ 0 −0.0176 0.0354 −0.0786 ⎟ ⎜ ⎟ 0.0071 ⎟ ⎜ 0 0.0364 0.0005 ×⎜ ⎟. ⎜ 0 −0.0131 −0.0029 −0.0007 ⎟ ⎜ ⎟ ⎜ 0 −0.0105 −0.0021 0.0006 ⎟ ⎝ ⎠ 0 0.0051 0.0001 0.0198 0 0.0027 −0.0008 −0.0010

0.0149 0.3779 −0.4755 −0.0718 0.1862 −1.2964 0.2267 −0.0411 −0.0562 0.0688 0.0962 0.0203 −0.1915 0.0012 0.1000 0.1287 0.0535 0.1106 2.0076 0.0726 −0.0611

−0.5229 0.2878 −0.0147 −0.0105 −0.0165 −0.0173 −0.0280

−0.0814 ⎟ ⎟ ⎟ ⎟ 0.4218 ⎟ ⎟ ⎟ ⎟ −0.5178 ⎟ ⎟ ⎟ ΛM (k)⎟ 0.0889 ⎟ ⎟ ⎟ ⎟ 0.0506 ⎟ ⎟ ⎠ ⎠ −0.1545 −0.1271

J.V. Salcedo, M. Martínez / Applied Numerical Mathematics 58 (2008) 1553–1571

1571

References [1] B. Bamieh, L. Giarré, Identification of linear parameter varying models, International Journal of Robust and Nonlinear Control 12 (2002) 841–853. [2] B. Bodenheimer, P. Bendotti, Optimal linear parameter-varying control design for a pressurized water reactor, in: Proceedings of the 34th Conference on Decision and Control, 1995, pp. 182–187. [3] B. Bodenheimer, P. Bendotti, M. Kanter, Linear parameter-varying control of a ducted fan engine, International Journal of Robust and Nonlinear Control 6 (9–10) (1996) 1023–1044. [4] S. Boyd, L. El Ghaoui, E. Feron, V. Balakrishnan, Linear Matrix Inequalities in System and Control Theory, Studies in Applied Mathematics, SIAM, Philadelphia, PA, 1994. [5] S. Dussy, L. El Ghaoui, Multiobjective bounded control of uncertain nonlinear systems: An inverted pendulum example, in: Control of Uncertain Systems with Bounded Inputs, Springer, Berlin, 1997, pp. 55–73. [6] L. El Gahoui, G. Scorletti, Control of rational systems using linear-fractional representations and linear matrix inequalities, Automatica 32 (9) (1996) 1273–1284. [7] L. El Ghaoui, S. Niculescu (Eds.), Advances in Linear Matrix Inequality Methods in Control, Advances in Design and Control, SIAM, Philadelphia, PA, 2000. [8] L. Giarre, D. Bauso, P. Falugi, B. Bamieh, LPV model identification for gain scheduling control: An application to rotating stall and surge control problem, Control Engineering Practice 14 (2006) 351–361. [9] L. Guzzella, A. Amstutz, Control of diesel engines, IEEE Control Systems Magazine 8 (9) (1998) 55–71. [10] L. Guzzella, C.H. Onder, Introduction to Modeling and Control of Internal Combustion Engine Systems, Springer, 2004. [11] L.H. Lee, K. Poolla, Identification of linear parameter-varying systems using nonlinear programming, Journal of Dynamical Systems, Measurement and Control 121 (1999) 71–78. [12] L. Ljung, System Identification. Theory for the User, second ed., Prentice-Hall, Englewood Cliffs, 1999. [13] M.C. Mazzaro, B.A. Movsichoff, R.S. Sánchez Peña, Robust identification of linear parameter varying systems, in: Proceedings of the American Control Conference, 1999, pp. 2282–2284. [14] M. Nemani, R. Ravikanth, B.A. Bamieh, Identification of linear parametrically varying systems, in: Proceedings of the 34th Conference on Decision and Control, 1995, pp. 2990–2995. [15] W.J. Rugh, J.S. Shamma, Research on gain scheduling, Automatica 36 (2000) 1401–1425. [16] J.V. Salcedo, X. Blasco, M. Martínez, J.V. García, Modelling and control in Simulink of a turbocharged diesel engine, in: XXII Jornadas de Automática (Spain), 2001 (in Spanish). [17] J.V. Salcedo, M. Martínez, J. Sanchis, X. Blasco, Design of GPC’s in state space, Automatika 42 (3–4) (2002) 159–167. [18] C.W. Scherer, LPV control and full block multipliers, Automatica 37 (2001) 361–375. [19] M. Sznaier, M.C. Mazzaro, An LMI approach to control-oriented identification and model (in)validation of LPV systems, IEEE Transactions on Automatic Control 48 (9) (2003) 1619–1624. [20] V. Verdult, M. Verhaegen, Subspace identification of multivariable linear parameter-varying systems, Automatica 38 (2002) 805–814. [21] M.G. Wassink, C. Van de Val, M. Scherer, O. Bosgra, LPV control for a water stage: Beyond the theoretical solution, Control Engineering Practice 13 (2005) 231–245.