Pascal Dufour
email: dufour@lagep.univ-lyon1.fr

Madiha Nadri
email: nadri@lagep.univ-lyon1.fr

Jun Qian
email: qian@lagep.univ-lyon1.fr

An explicit optimal input design for first order systems identification

Keywords: Optimal experiment design, input design, identification, linear systems, observers, predictive control

This paper focuses on the problem of closed loop online identification of the time constant in the single input single output (SISO) first order linear model. A new explicit approach for the simultaneous online optimal experiment design (OED) and model parameter identification is presented. Based on the observation theory and a model based predictive control (MPC) algorithm, this approach aims to solve an optimal control problem where input and output constraints may be specified. This constrained control objective aims to maximize the sensitivity of the model output with respect to the unknown model parameter (the time constant). The control law is derived explicitly offline and simple to be implemented: the input may be computed fast online while the unknown model time constant is estimated at the same time.

INTRODUCTION

When using dynamic models for simulation, control or optimization, all model parameters need to be numerically known. Among the dynamic models, the first order linear model is the simplest model: it is based on a static gain and a time constant. Methods to estimate the value of the parameters of a considered model concern the identification task (see Ljung (1999)), which relies on experimental data. OED may be useful if no data exist for identification, since it aims to design new experimental data (see Franceschini and Macchietto (2008)) which are often used offline for the identification procedure, hence decoupling these two procedures. This paper focuses on the coupled online OED and parameter estimation for such a first order linear model, where the time constant (that may vary during the time) has to be estimated. For such estimation problem, one may first think to apply a step input. But, this input in not sufficiently persistent (i.e., when the steady is reached in that case after some time, the output does not contain enough information on the dynamic, and the time constant can not be estimated). Pseudo-random binary sequence (PRBS) is a better open loop approach, since it features a rich input signal. However, open loop control does not allow to maintain the process output is a prescribed region, for example if the process has still to be used in the meantime for production. Therefore, an input design for identification based on closed loop control is required. In Shouche et al. (2002), the authors developed a constrained MPC based on an auto regressive model with external inputs (ARX) with a persistent input constraint, where the parameters are estimated online by regression. This led to a non-convex simultaneous online identification and control problem. The new method developed here for the SISO first order linear model is a particular development of a general approach (Qian et al. (2013(Qian et al. ( , 2014))) developed for multivariable nonlinear models. In these papers, an approach combining the closed loop online parameter identification with OED has been proposed for the input design. Based on a MPC, it aims to compute online the optimal input that maximize a norm of the sensitivities of the process outputs with respect to the unknown model parameters, which are estimated at the same time by an observer. The predictive controller handles input and output constraints while its cost function is based on the sensitivity criterion.

Here, based on this approach, for the SISO first order linear model, the problem is to design a controller and an estimator that both allow to estimate the unknown time constant (that may varies) and to maintain the process in a prescribed region. Based on the linearity of the model and assuming that the argument of the MPC is a single value (hence, the control horizon is tuned to one), an explicit controller can be derived offline, which avoid to solve any online optimization task. This paper is structured as follows: Section 2 deals with some recalls on the previous general approach of Qian et al. (2013Qian et al. ( , 2014) ) developed for multivariable nonlinear models. Section 3 contains the proposed algorithm for the constrained optimal closed loop online identification of the time constant in a SISO first order linear model. In section 4, the simulation results demonstrate the efficiency of this approach.

PRELIMINARIES

In Qian et al. (2013Qian et al. ( , 2014)) Online computations

• at each current time t k , get the measure of the process output y p (t k ) • then, based on the process input applied u(t k ) and y p (t k ), integrate the observer to get the current estimate of the unknown parameters • then, based on the measure y p (t k ), the parameter estimations, the nonlinear MPC constrained optimization problem is solved since, in a general case, and due to the nonlinearity of the problem, the search for the optimal control u max (t k ) to be applied at each sample time requires the online integration of the 2 prediction models used with the nonlinear MPC resolution. • then, u max (t k ) is applied from t k + to t k+1 , where is the time necessary for the whole loop (i.e, take the measure, integrate the observer, solve the constrained optimization problem). Hence, the main cost in time is related to the resolution of the constrained optimization problem involving model integrations, which might be important compared to the sampling time.

NEW ALGORITHM

Main ideas

For the previous approach, let us now consider the simplest SISO case: the well known first order model where x ∈ R is the measured state, u is the measured input, with one unknown parameter (the time constant T > 0), the known static gain G ∈ R , y p (t k ) is the measured process output at the current time t k considered as the initial time:

   ẋ(t) = -1 T x(t) + G T u(t), t > t k x(t = t k ) = x(t k ) = y p (t k ) (1)
In that case, the whole online numerical resolution may be simplified to get an explicit controller derived offline, hence reducing the online computational burden.

Offline synthesis The control law may be designed offline as follows:

• define the sensitivity dynamic model, in the time domain. Online computations

• at each current time t k , get the process measure y p (t k )

• then, based on the process input applied u(t k ) and y p (t k ), integrate the observer to get the current estimate of the time constant T (t k ) • then, based on the measure y p (t k ) and the estimation T (t k ), the optimal control to apply u(t k ) = u max (t k ) (constrained in magnitude) is computed, simply obtained from the evaluation of the explicit control law defined offline based on u min (t k ) and the input bounds. • then, u max (t k ) is applied from t k + to t k+1 , where is the time necessary for the whole loop (i.e, take the measure, integrate the observer, compute u max (t k )). Hence, the main cost in time is the observer integration: it is low in that case compared to the sampling time (if it is larger than the order of the ms). Therefore, the proposed algorithm for the SISO first order model is a very fast control algorithm.

Following the existing general method given in Qian et al. (2013,2014), the steps of the main ideas for this new input design are detailed in the following parts.

Sensitivity model

Using the considered SISO model ( 1), the sensitivity model state x s = ∂x ∂T is given by:

     ẋs (t) = 1 T 2 x(t) - 1 T x s (t) - G T 2 u(t), t > t k x s (t = t k ) = x s (t k ) x s (t = 0) = 0 (2)

Observer design

From the first order model (1), the augmented state is

x a (t) = [x(t) 1 T ] T in the augmented model: ẋa (t) = A a (y p (t), u(t))x a (t), t > t k ŷ(t) = C a x a (t), t > t k (3)
where:

A a (y p (t), u(t)) = 0 -y p (t) + Gu(t) 0 0 C a = [1 0] (4)
Remark 1. The system (3-4) is observable if (u(t), y p (t)) is a persistent input [START_REF] Besançon | Nonlinear observers and applications[END_REF]).

Different observers based on the augmented model (3-4) may be used to estimate the unknown model parameter T [START_REF] Besançon | Nonlinear observers and applications[END_REF]). A possible observer, which state is xa (t) = [x(t) 1 T (t) ] T , may be written here as an extended Kalman filter (EKF):

               ẋa (t) = A a (y p (t), u(t))x a (t) -ρ o S -1 o C T a (C a xa (t) -y p (t)), t > t k Ṡo (t) = -θ o S o (t) -A T a (t)S o (t) -S o (t)A a (t) + ρ o C T a C a , t > t k xa (t = t k ) = xa (t k ) = [y p (t k ) 1 T (t k ) ] S o (t = t k ) = S o (t k )
(5) where ρ o > 1, θ o > 0 and S o (t = 0) (a symmetric positive definite matrix) are the observer tuning parameters.

Sensitivity model step response

In the MPC framework, the optimization procedure aims to find the optimal u(t k ) to be applied during the next sample time through the maximization of the norm of the sensitivity dynamic model state. This requires here to evaluate x s (t) into the prediction horizon N p . Starting from the current time t = t k , let us apply a step u(t k ) for the process input in the MPC procedure. We use the Laplace transform formula, for a function f (t) starting from any initial condition f (t k ) at any initial time t k , is:

     L( df dt (t)) = sL(f (t)) -f (t k ), t > t k ⇐⇒ L( df dt (t)) = sF (s) -f (t k ), t > t k (6)
Therefore, applying the Laplace transform to the system (1-2), at the current time t k leads to:

X s (s) = x s (t k ) s + 1 T + (y p (t k ) -GsU (s)) T 2 (s + 1 T ) 2 , t > t k . ( 7
)
where X s (s) is the Laplace transform of x s (t). Based on the Laplace transform of the input step function fed into (7), the response of the sensitivity model over the prediction horizon N p is (in the Laplace domain):

X s (s) = x s (t k ) s + 1 T + (y p (t k ) -Gu(t k )) T 2 (s + 1 T ) 2 , t > t k (8)
Hence, the predicted sensitivity response into the future requires the current process measure y p (t k ), the unknown model time constant T and of course the input value u(t k ).

Based on the usual inverse Laplace transform formula, the step response x s (t) into the prediction horizon N p may be written in the time domain, where the unknown parameter T is replaced by the last current estimation from the observer at t k ( T (t k )):

x s (t) = ...

x s (t k ) + (y p (t k ) -Gu(t k ))(t -t k ) T (t k ) 2 e -(t-t k ) T (t k ) , t > t k (9)
From ( 9), let us notice two well known particular cases for the first order system which have to be avoided to be able to estimate the time constant:

(1) In open loop, starting from any initial condition (y p (0), x s (0) = 0), if u(t) is kept constant, after a certain time, the sensitivity tends to 0. (2) Moreover, in the previous case, if the initial condition is a steady state (hence y p (t) -Gu(t) = 0, ∀t ≥ 0 and since x s (0) = 0), then the sensitivity is 0 ∀ t .

Constrained optimal input design

In the present case, the MPC formulation of Qian et al. (2013,2014) is written as:

   max u inf ≤u(t k )≤u sup J(u(t k )) = t k+Np t k (x s (t)) 2 dt where x s (t) is obtained from (9) (10)
where the control bounds u inf and u sup are adjusted offline according to the prescribed constraints on u and y p :

u inf,u ≤ u ≤ u sup,u , ∀t > 0 y inf ≤ y p ≤ y sup , ∀t > 0 (11)
In the following, if there is no such input constraints, u inf,u is set to -∞ and u sup,u is set to +∞. If there is no such output constraints, y inf is set to -∞ and y sup is set to +∞. In there is neither such input constraints nor output constraints, the problem is not well defined. Since any step response of a first order system is monotonous and takes its maximum (or minimum) value at the steady state, and taking account of the sign of the static gain (it may be positive or negative) then the control bounds u inf and u sup in (10) may be simply tuned offline from ( 11) and the known static gain G:

             u sup,y = max( y inf G , y sup G ) u sup = min(u sup,y , u sup,u ) u inf,y = min( y inf G , y sup G ) u inf = max(u inf,y , u inf,u ) (12)
Due to the input to output linearity of the sensitivity model, the sensitivity involved in the integration in ( 10) is quadratic in its optimization argument u(t k ):

                                             max u inf ≤u(t k )≤u sup J(u(t k )) = t k+Np t k f 2 (t)u(t k ) 2 + f 1 (t)u(t k ) + f 0 (t)dt
where:

f 2 (t) = e -2(t-t k ) T (t k ) G 2 T (t k ) 4 (t -t k ) 2 ∈ R + f 1 (t) = e -(t-t k ) T (t k ) ... -2(x s (t k ) + y p (t k ) T (t k ) 2 (t -t k )) G T (t k ) 2 (t -t k )e -(t-t k ) T (t k ) f 0 (t) = e -(t-t k ) T (t k ) x s (t k ) + y p (t k ) T (t k ) 2 (t -t k )e -(t-t k ) T (t k ) (13)
Due to the positive sign of f 2 (t), the problem ( 13) is convex and may be rewritten showing characteristics of this parabola:

max u inf ≤u(t k )≤u sup J(u(t k )) = ... ...c ∞ (t k )(u(t k ) -u min (t k )) 2 + J min (t k ) (14)
where c ∞ (t k ) is a positive real (obtained from the integration of f 2 (t) in ( 13)), and where u min (t k ) is the worst control u(t k ) to apply, since it minimizes at t k the criteria J(.) (where J(.) = J min (t k ) ≥ 0, since J is the L 2 norm of x s ). Hence, u min (t k ) has to be determined through the first order optimality condition:

∂J ∂u(t k ) = 0 at u(t k ) = u min (t k ) (15)
which leads to:

u min (t k ) = - t k+Np t k f 1 (t)dt 2 t k+Np t k f 2 (t)dt (16)
This leads to determine u min (t k ) by the integration of the two dynamic models, starting at each t k from the two initial conditions y p (t k ) and x s (t k ) and with the estimation T (t k ):

           u min (t k ) = 1 G y p (t k )... ... + T (t k ) 2 G t k+Np t k (t -t k )e -2(t-t k ) T (t k ) dt t k+Np t k (t -t k ) 2 e -2(t-t k ) T (t k ) dt x s (t k ) (17)
By pursuing the integration and replacing t k+Np -t k by the product N p T s (where T s is the constant sampling time), this leads to:

                       u min (t k ) = 1 G T (t k ) 2 γ(.)x s (t k ) + y p (t k ) with γ(N p , T s , T (t k ), ν(N p , T s , T (t k ))) = 2N p T s ν(.) + T (t k )(ν(.) -1) 2(N p T s ) 2 ν(.) + 2 T (t k )N p T s ν(.) + T (t k ) 2 (ν(.) -1) and ν(N p , T s , T (t k )) = e -2Np Ts T (t k ) (18)
Due to the symmetry property of such a parabola with respect to the vertical line, which is passing through (u min (t k ),J min (t k )), one may conclude with the optimal control law for u max (t k ):

u max (t k ) =            u inf if u min (t k ) > u inf + u sup -u inf 2 u sup
otherwise where: u sup and u inf are defined in ( 12) and u min (t k ) is defined in (18) (19) Therefore, the optimal control law u max (t k ) may be explicitly defined offline and is a bang bang control which value is evaluated online with the process measure and the time constant estimation, and is either the minimum of the maximum input value allowed. In as to be noted that the proposed solution is derived using the first order condition, and the constrained are enforced later resulting in bangbang control. In a general framework, this approach is not necessarily optimal. But here, since the first order condition condition is used to find the worst control (and not the best one), the optimal control is really found (and not a suboptimal solution).

SIMULATION RESULTS

Simulations allow to compare 2 controller performances:

• case OL: an open loop PRBS input, • case CL: the optimal closed loop controller with input constraints and output constraints based on the algorithm (19).

The time constant T p of the simulated process is assumed to vary in 9 periods, between 3 values (50, 80 and 100s) as step functions, with 3 different length for each step (1000, 500 and 250s). Runs are done with Matlab1 under the following conditions (table 1):

T (0) 20 (s) Controller [u inf,u , u sup,u ] [0, 1] [y inf , y sup ] [0, 4]
After some trial-errors, the prediction horizon N p is tuned at each t k according to:

N p (t k )T s = T (t k ) Simulations results are summarized in table 2: they allow to see that for the 2 cases, in spite of a 300 % initial error given by the initial condition of the time constant fed into the observer, the estimation of the time varying process time constant is done with a good accuracy for each step. But the estimations is better in the closed , T f inal is 3.6%, see Fig. 1) rather than 4.62% in open loop where it is more noisy (Fig. 2). Let us notice that in open loop, at the beginning, a very large estimation error occurs, which is not the case in closed loop (the short period where it occurs is not in T f inal 100 , T f inal ). The drawback of the OL is that the output is uncontrolled (ranging between 0 and 6.14 (Fig. 3)). If such output values are not acceptable, it is preferable to use the proposed closed loop control algorithm CL with both input constraints and output constraints (Fig. 4) where the output ranges between 0.05 and 3.65 (hence always within the bounds of the output constraints). Such output constraints impact u inf and u sup and hence the control values u min (t) and u max (t) applied (see Fig. 5). The sensitivity maximized in CL (see Fig. 6) is 3 times higher in average than in OL (see Fig. 7). In this paper, an explicit bang bang controller has been proposed for the estimation of the time constant in a first order linear model. This input design relies on an observer (solved online with the input and output measures) and a constrained MPC problem, which exact optimal solution is derived explicitly offline. Hence, online computational cost is very low since it concerns mostly the observer integration. Perspectives are concerned with the optimal tuning of the control horizon and the robustness of this approach with respect to output noise and model uncertainties.

Fig. 1 .

 1 Fig. 1. Time constant: time varying target and estimation (with closed loop control with constraints on u and y)

Fig. 3 .

 3 Fig. 3. Output (with open loop control with constraints on u)

Fig. 6 .

 6 Fig. 6. Parameter sensitivity (with closed loop control with constraints on u and y)

  • define an augmented dynamic model, where the state contains the model state and the unknown time constant. It is used to define an observer that aims to estimate the unknown time constant. • get the predicted step response of the sensitivity model into the prediction horizon N p . It is based on the Laplace transform applied to the process model and to the sensitivity model: hence, the formula is parametrized for any initial conditions of the model state x(t k ) and the sensitivity model state x s (t k ) • due to the convexity of the constrained quadratic problem, get the worst control law u min (t

k ) (that minimize the MPC cost function) and therefore the persistent control law u max (t k ) (that maximize the MPC cost function as requested by the method) to be applied at each time t k .

Table 1 .

 1 Simulation conditions

		Parameter	Value (unit)
	Time	T f inal	5250 (s)
		Ts	1 (s)
	Model and Process	G	10 (-)
		yp(0)	0 (-)
		Tp(0)	80 (s)
	Observer	ρo	1.02 (-)
		θo	0.05 (-)
		So	0.01 0 0 0.01	(-)

Table 2 .

 2 Results analysis

		Criteria	Input
			OL	CL
	Maximum value taken by yp	6.14	3.655
	Minimum value taken by yp	0.0	0.050
	Mean sensitivity over
		T f inal 100 , T f inal ( 1 s )	0.004	0.011
	Mean time constant estimation error
	over	T f inal 100 , T f inal (%)	4.62	3.60

www.mathworks.com

company Acsystème which is gratefully acknowledged for the funding. The french ministry of higher education and research, for the financial support of this CIFRE PhD thesis 2011/0876,


Chouaib Afri
email: afri@lagep.univ-lyon1.fr

Vincent Andrieu
email: vincent.andrieu@gmail.com

Laurent Bako
email: laurent.bako@ec-lyon.fr

Pascal Dufour
email: dufour@lagep.univ-lyon1.fr

Identification of linear systems with nonlinear Luenberger Observers

The design of a nonlinear Luenberger observer for a parametrized linear system is studied. From an observability assumption of the system, the existence of such an observer is concluded. In a second step, a constructive novel algorithm for the identification of multi-input multi-output linear systems is suggested and is implemented on a second order system.

I. INTRODUCTION

In this note, a parametrized linear system described by the following equation is considered:

ẋ = A(θ )x + B(θ )u , y = C(θ )x, (1)
where θ in Θ ⊂ R q is a vector of unknown parameters and Θ is a known set, u in R m is a controlled input (in L ∞ loc (R + )). The state x is in R n and y is the measured outputs in R p . Mappings A : Θ → R n×n , B : Θ → R n×m and C : Θ → R p×n are known C 1 matrix valued functions.

In the following, an estimation problem is considered: The design of an observer to estimate the state and the unknown parameters of the system from the knowledge of y is seek. In other word, an asymptotic observer for the extended (nonlinear) n + q dimensional system ẋ = A(θ )x + B(θ )u , θ = 0 , y = C(θ )x

(2)

has to be designed. The nonlinear Luenberger methodology introduced in [START_REF] Luenberger | Observing the state of a linear system[END_REF] (see also [START_REF] Shoshitaishvili | Singularities for projections of integral manifolds with applications to control and observation problems[END_REF], [START_REF] Kazantzis | Nonlinear observer design using Lyapunov's auxiliary theorem[END_REF], [START_REF] Kreisselmeier | Nonlinear observers for autonomous Lipschitz continuous systems[END_REF], [START_REF] Andrieu | On the existence of Kazantzis-Kravaris / Luenberger Observers[END_REF], [START_REF] Andrieu | Convergence speed of nonlinear luenberger observers[END_REF]) is a method which permits to design an observer based on weak observability assumptions. Following this approach, the first step is to design a C 1 function (x, θ , w) → T (x, θ , w) such that the following partial differential equation (PDE) is satisfied:

∂ T ∂ x (x, θ , w)[A(θ )x + B(θ )u] + ∂ T ∂ w (x, θ , w)g(w, u) = ΛT (x, θ , w) + LC(x, θ )x , (3)
where Λ is a Hurwitz squared matrix, L a matrix and g is a controlled vector field. Dimensions of the matrices and of the vector field g must be taken consistently. This will be precisely defined in the following. The interest on this mapping is highlighted if (z(•), w(•)), the solution initiated from (z 0 , w 0 ) of the dynamical system, is considered: ż = Λz + Ly , ẇ = g(w, u) .

Indeed, assuming completeness (of the w part of the solution), for all positive time t: ˙ z(t) -T (x(t), θ , w(t)) = Λ(z(t) -T (x(t), θ , w(t))) .

Hence, due to the fact that Λ is Hurwitz, *asymptotically get lim t→+∞ |z(t) -T (x(t), θ , w(t))| = 0 .

In other words, z provides an estimate of the function T . The second step of the Luenberger design is to left invert the function T in order to reconstruct the extended state (x, θ ) from the estimate of T . Hence, a mapping T * has to be constructed such that

T * (T (x, θ , w), w) = (x, θ ) . (5)
Of course, this property requires the mapping T to be injective. The final observer is simply

ż = Λz + Ly , ẇ = g(w, u) , ( x, θ ) = T * (z, w) .
A particularly interesting feature of this observer is that its convergence rate can be made as large as requested (see [START_REF] Andrieu | Convergence speed of nonlinear luenberger observers[END_REF]). In this paper, this strategy is adopted to suggest a solution to the state and parameters estimation for system [START_REF] Andrieu | On the existence of Kazantzis-Kravaris / Luenberger Observers[END_REF]. This paper can be seen as an extension of the result of [START_REF] Praly | A new observer for an unknown harmonic oscillator[END_REF] in which a nonlinear Luenberger observer is constructed for a harmonic oscillator which fits in the class of system studied.

As will be seen in Section IV, our approach can be employed to identify the state and all parameters describing an unknown linear model. This topic has been widely studied in the literature (see the books [START_REF] Ioannou | Robust adaptive control[END_REF], [START_REF] Marino | Nonlinear control design: geometric, adaptive and robust[END_REF], [START_REF] Ks Narendra | Stable Adaptive Systems[END_REF]). Adaptive identifiers can be traced back to G. Kreisselmeier in [START_REF] Kreisselmeier | Adaptive observers with exponential rate of convergence[END_REF]. In this paper three algorithms were given for identification of a SISO system in a specific controllable and observable representation. This work has then been extended in many directions to allow time varying matrices and MIMO systems (see for instance [START_REF] Bastin | Stable adaptive observers for nonlinear time-varying systems[END_REF], [START_REF] Marino | Adaptive observers with arbitrary exponential rate of convergence for nonlinear systems[END_REF], [START_REF] Zhang | Adaptive observer for multiple-input multiple-output (mimo) linear time-varying systems[END_REF]).

The algorithm given in Section IV gives a new approach to address the same problem. It has the advantage to allow for a prescribed convergence rate for MIMO systems. Moreover, its structure is relatively simple (see ( 27)) and its dimension significantly smaller than that of available algorithms. One of the main difference is that, in contrast to all other available approaches, no optimization step is required.

The paper is organized as follows: In Section II, the existence of a mapping T which solves the partial differential equation ( 3) is discussed. Section III is devoted to the study of the injectivity of the mapping T assuming some observability properties. In Section IV, a left inverse of the mapping T is constructed to get the observer when considering a specific structure for the matrices A, B and C. This leads to a novel solution for the identification of MIMO linear time invariant systems.

Notations:

• Given a matrix A in R n×n , σ (A) denotes its spectrum.

• 1 n denotes the n dimensional real vector composed of 1.

• I n denotes the n dimensional identity matrix.

• ⊗ stands for the Kronecker product between two matrices.

• Given a C j function u: ū( j) (t) = u . . . u ( j) .

• For a vector or a matrix | • | denotes the usual 2-norm.

• Given a function ρ : R + × T → R + the notation ρ(s,t) = o(s) is used if for all ε > 0, there exists s 0 > 0 such that for all (s,t) in (0, s 0 ) × T such that:

ρ(s,t) ≤ εs .

• Given a set C, Cl(C) is its closure.

II. EXISTENCE OF THE MAPPING T

In [START_REF] Andrieu | On the existence of Kazantzis-Kravaris / Luenberger Observers[END_REF], it is shown that, in the autonomous case the existence of the mapping T , solution of the PDE (3), is obtained for almost all Hurwitz matrices Λ. In the context of the controlled system (2), the same type of result still holds. Moreover, an explicit solution of the PDE (3) may be given.

Theorem 1 (Existence of T ): Let r be a positive integer. For all r-uplet of negative real numbers (λ 1 , . . . , λ r ) such that

λ i / ∈ θ ∈Θ σ {A(θ )} , i = 1, . . . , r, (6)
there exists a linear in x function T : R n × Θ × R mr → R pr solution to the PDE (3) with Λ = Diag{λ 1 , . . . , λ r }⊗I p , L = 1 r ⊗ I p , and g : R rm × R m → R rm defined as g(w, u) = Λ m w + L m u with Λ m = Diag{λ 1 , . . . , λ r } ⊗ I m , L m = 1 r ⊗ I m . Proof: Keeping in mind that the spectrum of Λ and A(θ ) are disjoint (see [START_REF] Kreisselmeier | Adaptive observers with exponential rate of convergence[END_REF]), the matrix M i (θ ) in R p×n can be defined for all i in {1, . . . , r}:

M i (θ ) = C(θ )(A(θ ) -λ i I n ) -1 .
For all i = 1, . . . , r, let T i : R n × Θ × R m → R p and the vector field g i : R m × R m → R m be defined as:

T i (x, θ , w i ) = M i (θ )[x -B(θ )w i ] , g i (w i , u) = λ i w i + u , (7)
It can be noticed that T i is solution to the PDE

∂ T i ∂ x (x, θ , w i )[A(θ )x + B(θ )u] + ∂ T i ∂ w i (x, θ , w i )g i (w i , u) = λ i T i (x, θ , w i ) +C(θ )x .
Hence, the solution of the PDE (3) is simply taken as

T (x, θ , w) = T 1 (x, θ , w 1 ) . . . T r (x, θ , w r ) . (8)
This ends the proof.

Remark 1: Note that in the particular case in which the system is autonomous and there is only one output, the mapping T is given as

T o (x, θ ) = M(θ )x , M(θ ) =    M 1 (θ ) . . . M r (θ )    . ( 9
)
This matrix is solution to the following parametrized Sylvester equation

M(θ )A(θ ) = ΛM(θ ) + LC(θ ) . (10)
Hence, taking r = n, the well known Luenberger observer introduced in [START_REF] Luenberger | Observing the state of a linear system[END_REF] in the case of autonomous systems is recovered. Note however that, here, the injectivity is more involved than in the context of [START_REF] Luenberger | Observing the state of a linear system[END_REF] since θ is unknown.

Remark 2: Note that if the set Θ is bounded, then it is ensured that there exists (λ i )'s which satisfy equation [START_REF] Kreisselmeier | Adaptive observers with exponential rate of convergence[END_REF]. Indeed, if Θ is bounded, then the set ( θ ∈Θ σ {A(θ )}) is a bounded set. This can be obtained from the fact that each eigenvalues λ in σ {A(θ )} is a zero of the characteristic polynomial:

λ n + µ 1 (θ )λ n-1 + • • • + µ n-1 (θ )λ + µ n (θ ) = 0, (11)
where µ i (θ ) are continuous function of θ . Boundedness of Θ together with the continuity of the µ i 's imply that there

is c > 0 such that |µ i (θ )| ≤ c ∀i ∈ {1, . . . , n}, ∀θ ∈ Θ. As a consequence if |λ | > 1, we must have |λ | ≤ n ∑ j=1 |µ j (θ )||λ | 1-j ≤ c 1 -1/|λ | which hence implies that |λ | ≤ c + 1.

III. INJECTIVITY

As seen in the previous section, it is known that if the following dynamical extension is considered:

ż = Λz + Ly , ẇ = Λ m w + L m u (12)
with z in R rp and w in R rm , then it yields that along the solution of the system (2)-( 12), equation ( 4) is true. Consequently, T (x, θ , w) defined in ( 7)-( 8) is asymptotically estimated. The question that arises is whether this information is sufficient to get the knowledge of x and θ . This is related to the injectivity property of this mapping. As shown in [START_REF] Andrieu | On the existence of Kazantzis-Kravaris / Luenberger Observers[END_REF], in the autonomous case this property is related to the observability of the extended system [START_REF] Andrieu | On the existence of Kazantzis-Kravaris / Luenberger Observers[END_REF]. With observability, it is sufficient to take r large enough to get injectivity. Here, the same type of result holds if it is assumed an observability uniform with respect to input in a specific set.

In order to simplify the presentation, the SISO case is considered. The following strong observability assumption is made:

Assumption 1 (Uniform differential injectivity): There exist two bounded open subsets C θ ⊂ Θ, C x ⊂ R n , an integer r and U r a bounded subset of R r such that the mapping

H r (x, θ , v) = H r (θ )x + r ∑ j=1 S j H r (θ )B(θ )v j-1 , with H r (θ ) =      C(θ ) C(θ )A(θ ) . . . C(θ )A(θ ) r-1      v = (v 0 , . . . , v r-1
) and S is the shift matrix operator such that for all s = (s 1 , . . . , s r ), S × s = (0, s 1 , . . . s r-1 ) is injective in (x, θ ), uniformly in (x, θ , v) ∈ Cl(C θ )×Cl(C x )×U r and full rank. More precisely, there exists a positive real number L such that for all (x, θ ) and

(x * , θ * ) in Cl(C θ ) × Cl(C x ) and all v in U r |H r (x * , θ * , v) -H r (x, θ , v)| ≥ L (|x -x * | + |θ -θ * |) .
The following result establishes an injectivity property for large eigenvalues of the observer.

Theorem 2: Assume Assumption 1 holds. Let u(•) be a bounded C r-1 ([0, +∞]) function with bounded r -1 first derivatives, i.e. there exists a positive real number u such that

| ū(r-1) (t)| ≤ u , ∀t ≥ 0 . ( 13
)
For all r-uplet of distinct negative real numbers ( λ1 , . . . , λr ), for all positive time τ and for all w 0 in R r , there exist two positive real numbers k * and L T such that for all k > k * the mapping defined in ( 7)-( 8) with λ i = k λi , i = 1, . . . , r satisfies the following injectivity property in

C = C x × C θ . For all t 1 ≥ τ, if ū(t 1
) is in U r then for all (x, θ ) and (x * , θ * ) in C x × C θ the following inequality holds:

|T (x, θ , w(t 1 )) -T (x * , θ * , w(t 1 ))| ≥ L T k r (|x -x * | + |θ -θ * |) (14)
where w(•) is the solution of the w dynamics in (12) initiated from w 0 .

Proof: First of all, picking k sufficiently large implies that the matrix M i which satisfies

M i (θ ) = 1 k λi C(θ ) 1 k λi A(θ ) -I n -1
is well defined. On another hand, assume that k is sufficiently large such that for all i in [1, r],

k| λi | ≥ |A(θ )| , ∀θ ∈ C θ . (15)
This implies that, for θ in C θ :

M i (θ ) = - r ∑ j=1 1 (k λi ) j C(θ )A(θ ) j-1 + R i (θ ) , with R i (θ ) = - +∞ ∑ j=r+1 C(θ )A(θ ) j-1 (k λi ) j .
Note moreover that R i is a C 1 function which satisfies for θ in C θ :

|R i (θ )| = o 1 k r , ∂ R i ∂ θ (θ ) = o 1 k r . (16)
Moreover, M i can be stated in the form which will be useful in the following: Let K be the matrix in R r×r defined as

K = Diag 1 k , . . . , 1 k r and Ṽi = 1 λi 1 λ 2 i . . . 1 λ r i
Note that:

M i (θ ) = -Ṽi K H r (θ ) + o 1 k r .
On another hand, for all t, since u is C r and w(•) being solution of (12) (in the SISO case), one get:

w ( j) i (t) = (k λi )w ( j-1) i (t) + u ( j-1) (t) , j = 1, . . . , r .
which implies that w i satisfies:

w i (t) = - r-2 ∑ j=0 u ( j) (t) (k λi ) j+1 + R wi (t)
where: R wi (t) = w

(r-1) i (t) (k λi ) r-1 = exp(k λi t)w (r-1) i (0) (k λi ) r-1 + t 0 exp((k λi )(t -s)) u (r-1) (k λi ) r-1 ds
and with:

w (r-1) i (0) = (k λi ) r-1 w i (0) + r-1 ∑ j=1 (k λi ) r-1-j u ( j) (0) .
Hence, with [START_REF] Shoshitaishvili | Singularities for projections of integral manifolds with applications to control and observation problems[END_REF], it yields that for all t :

|R wi (t)| ≤ C exp(k λi t) + u (k λi ) r
where C is a positive real number which depends on w i (0) and (u(0), . . . , u (r-2) (0)). Keeping in mind that λi is negative, when t is larger than τ > 0, the previous inequality satisfies

|R wi (t)| ≤ R wi0 (k) = o 1 k r , ∀t ≥ τ ,
where R wi0 depends on k but not on t. This implies that by collecting terms of higher order in 1 k in a function denoted R MBi :

M i (θ )B(θ )w i (t) = Ṽi K r ∑ j=1 S j H r (θ )B(θ )u ( j-1) (t)+R MBi (θ ,t)
and with (16) and using the fact that C θ and u ( j) (t) are bounded, for all t ≥ τ:

1 |R MBi (θ ,t)| = o 1 k r , ∂ R MBi ∂ θ (θ ,t) = o 1 k r .
Finally:

T i (x, θ , w i (t)) = Ṽi K H r (x, θ , ū(r-1) (t)) + R Ti (x, θ ,t)
and for all (x, θ ) in C x × C θ which is bounded and all t ≥ τ, it yields:

∂ R Ti (x, θ ,t) ∂ x = o 1 k r , ∂ R Ti (x, θ ,t) ∂ θ = o 1 k r .
This implies:

T (x, θ , w i (t)) = Ṽ K H r (x, θ , ū(r-1) (t))+R T (x, θ ,t) , (17)
where Ṽ in R r×r is the Vandermonde matrix defined as:

Ṽ =    1 λ1 •• 1 λ r 1 : : 1 λr •• 1 λ r r    ,
1 In this equation and in the following the notation o

1 k r is uniform in t ≥ τ.
and R T is a C 1 function which satisfies for all (x, θ

) in C x × C θ and t ≥ τ: ∂ R T (x, θ ,t) ∂ x = o 1 k r , ∂ R T (x, θ ,t) ∂ θ = o 1 k r .
Hence the mapping R T is globally Lipschitz with a Lipschitz constant in o 1 k r . Hence, it is possible to find k 0 such that for all k ≥ k 0 and all quadruples (x, x * , θ , θ * ) in C 2

x × C 2 θ , for all t ≥ τ:

|R T (x, θ ,t) -R T (x * , θ * ,t)| ≤ L 2 Ṽ -1 k r (|x -x * | + |θ -θ * |) . (18)
It can be shown that the result holds with this value of k 0 . Employing (17), it yields that, for all t:

|T (x, θ , w(t)) -T (x * , θ * , w(t))| ≥ Ṽ K H r (x, θ , ū(r-1) (t)) -H r (x * , θ * , ū(r-1) (t)) -|R T (x, θ ,t) -R T (x * , θ * ,t)| . ≥ H r (x, θ , ū(r-1) (t)) -H r (x * , θ * , ū(r-1) (t)) | Ṽ -1 ||K -1 | -|R T (x, θ ,t) -R T (x * , θ * ,t)| .
Consider now t 1 ≥ τ, the last term of the previous inequality can be lower-bounded by (18). Moreover, if u(t 1 ) . . . , u (r-1) (t 1 ) is in U r , the other term can be lowerbounded based on Assumption 1 and the result follows.

IV. CASE OF IDENTIFICATION PROBLEMS

A. State and parameter estimation algorithm

A concrete algorithm for computing simultaneously the state and the parameter vector θ of the linear system ( 1) is now derived. Such estimates can be obtained as:

( x(t), θ (t)) = arg min x,θ |T (x, θ , w(t)) -z(t)| 2 (19)
where z(t) is the value obtained by integrating the dynamical system [START_REF] Praly | A new observer for an unknown harmonic oscillator[END_REF] starting from an arbitrary user-defined initial value (z(0), w(0)). The formulation (19) can be seen as one way of inverting the mapping T (see ( 5)). A distinctive feature of the optimization problem (19) which is worth pointing out is that all data involved in it are signals at time t.

Hence there is no explicit use of the data relevant to the past. Solving the problem (19) can be very challenging, because the mapping T is in general nonlinear and can be highly complex. As a consequence, the optimization problem ( 19) is non-convex so that solving it directly might require an iterative searching scheme. One simple way to tackle it is as follows. Let im T ⊂ R n z , n z = pr, denote the range of the map T . Assume that z(t) ∈ im T for some t. Such an assumption seems quite reasonable since z(t) converges exponentially to im T as stated by [START_REF] Ioannou | Robust adaptive control[END_REF]. Under this assumption, the estimates x(t) and θ (t) can be searched such that:

T ( x(t), θ (t), w(t)) = z(t)
But explicit determination of the expression of those estimates is still missing. For this purpose, let us consider particular canonical structures for the matrix-valued functions A, B, C in (1) :

A(θ ) = A(θ a ) =        -θ a1 I p I p 0 •• 0 -θ a2 I p 0 I p •• 0 : : : . . . : : : : I p -θ ad I p 0 0 •• 0        ∈ R n×n , C(θ ) = C = I p 0 • • • 0 ∈ R p×n , B(θ ) = B(θ b ) = θ b1 • • • θ bm ∈ R n×m , ( 20
) where n = d p, θ bi ∈ R n , i = 1, . . . , m, and θ = [θ a , θ b ] is in R d+mn with θ a = θ a1 •• θ ad ∈ R d , θ b = θ b1 •• θ bm ∈ R mn .
Note that assuming such structures (20) for A, B,C is without loss of generality: any input-output behavior of a linear MIMO system can be described with a model of this structure. Such a realization is however not necessarily minimal but it is observable for any vector θ . The problem of identifying canonical representations of linear systems has been deeply investigated in the last three decades employing mainly standard identification techniques (see for instance [START_REF] Ioannou | Robust adaptive control[END_REF] and references therein). Our approach offers a new algorithm to address it. The interest of this one comes from the fact that the dimension of the proposed observer may not be too high and that the convergence rate is tunable (see [START_REF] Andrieu | Convergence speed of nonlinear luenberger observers[END_REF]). Let the mapping T be as in ( 7)- [START_REF] Luenberger | Observing the state of a linear system[END_REF].

z i (t) = M i ( θ ) x(t) -B( θ )w i (t) (21)
To carry out the calculations, it is convenient to observe that the matrices A(θ a ) and C in (20) can be written as:

A(θ a ) = A c ⊗ I p , C = e T 1 ⊗ I p ,
where ⊗ refers to the Kronecker product,

e 1 = [1 0 • • • 0] T ∈ R d and A c =      -θ a1 1 • • • 0 . . . . . . . . . . . . -θ a(d-1) 0 • • • 1 -θ ad 0 • • • 0      ∈ R d×d .
Exploiting these expressions and doing a little manipulation through Kronecker algebra lead to

M i (θ ) = (e T 1 ⊗ I p ) (A c -λ i I d ) -1 ⊗ I p . (22)
Let

J i =      -λ i 1 • • • 0 . . . . . . . . . . . . 0 • • • -λ i 1 0 • • • 0 -λ i      ∈ R d×d , then: A c -λ i I d = J i -θ a e T 1 .
Applying the Sherman-Morrison-Woodbury formula, one gets:

(A c -λ i I d ) -1 = J -1 i + J -1 i θ a e T 1 J -1 i 1 -e T 1 J -1 i θ a ( 23
)
where it is assumed that 1e T 1 J -1 i θ a = 0. Indeed it can be easily checked that this is true provided the assumption (6) holds. Combining ( 22) and (23) gives us:

(1e T 1 J -1 i θ a )M i (θ ) = (e T 1 J -1 i ) ⊗ I p , which, together with (21), reveals that:

(1 -e T 1 J -1 i θa (t))z i (t) = (e T 1 J -1 i ⊗ I p ) x(t) -(e T 1 J -1 i ⊗ I p )B( θb )w i . Setting V i = 1 λ i . . . 1 λ d i ∈ R d , it can be verified that V T i = e T 1 J -1 i .
Rearranging the previous expression, one gets:

z i = z i V T i V T i ⊗ I p -(w T i ⊗V T i ) ⊗ I p P i (z i ,w i )   θa x θb   (24)
Finally, the following linear equation for ( θa , θb , x) must be solved:

z(t) = P(z(t), w(t))   θa (t) x(t) θb (t)   (25)
where P(z, w) = P 1 (z, w)

• • • P r (z, w) ∈ R pr×(n+q) with q = nm + d. A closed-form expression of the estimate is then:   θa (t) x(t) θb (t)
  = (P(z, w) P(z, w)) -1 P(z, w) z(t) (26) whenever the matrix P is full rank. However, the solution (26) can be efficiently computed only if the matrix P(z(t), w(t)) has full column rank at any time t . This is indeed related to the injectivity of the map T and more precisely to its rank. Indeed, let us give the following proposition: Proposition 1: Assume that the matrix-valued functions A, B, C in system (1) have the canonical structure (20). Assume moreover that condition (6) holds. Then, for all (z, x, θ , w) such that z = T (x, θ , w) and such that:

rank ∂ T ∂ (x, θ ) (x, θ , w) = d + m(n + 1)
then the matrix P defined in (24)-(25) has full column rank.

Proof: The proof follows from equation (25) in the case in which z i = T i (x, θ , w). Indeed, differentiating with respect to (x, θ ) this equation yields to:

∂ T i ∂ (x, θ ) (x, θ , w) = ∂ T i ∂ (x, θ ) (x, θ , w)V i θ a +P i (T i (x, θ , w), w)
This implies that

[1 -V i θ a ] ∂ T i ∂ (x, θ ) (x, θ , w) = P i (T i (x, θ , w), w)
Hence:

Diag(1 -V T 1 θ a , ..., 1 -V T r θ a ) ∂ T ∂ (x, θ ) (x, θ , w) = P(z, w)
Again, since condition (6) holds, it yields that Diag(1 -V T 1 θ a , ..., 1 -V T r θ a ) is invertible and consequently P is full rank if ∂ T ∂ (x,θ ) (x, θ , w) is. Consequently, since it is known that z converges asymptotically to im T , if T is injective and full rank, the observer given in (26) is well defined after a certain time. However, in the transient period, there is no guarantee that this matrix is well defined. This problem has been deeply investigated in [START_REF] Praly | A new observer for an unknown harmonic oscillator[END_REF] considering an autonomous second order system with only one parameter.

B. Numerical illustration

In order to show the performances of the observer ( 12)-(26), a second order linear system in the form (1) is simulated in Matlab where:

A(θ ) = θ 1 1 θ 2 0 ; B(θ ) = θ 3 θ 4 ; C(θ ) = ( 1
0 ) θ i 's, the initial value x(0) and the input signal are given in Table I.

θ 1 = θ 2 = -1; θ 3 = θ 2 = 0.5 x(0) = (x 1 (0) x 2 (0)) = (1 -1)
u(t) = sin(t) + sin(2.2t) + sin(4.5t) + sin(6.7t) + sin(7t) It is assumed that the state (x 1 , x 2 ) and the parameters (θ 1 , θ 2 , θ 3 , θ 4 ) are unknown and are estimated with:

       ż(t) = Λz(t) + Ly(t) ẇ(t) = Λw(t) + Lu(t) ( x(t), θ (t)) = P T (w(t), z(t))P(w(t), z(t)) -1 × P T (w(t), z(t))z(t) (27)
where, P is obtained from (24). Λ, L and the initial configuration of z(0) and w(0) are given in Table II. The ordinary differential equations in (27) are solved with the Matlab lsim function. As is proved in Theorem 2 and Proposition 1, the rank of the matrix P is related to observability property and to the number of λ i 's considered. In this simulation, there are r = 7 eigenvalues which is more than the least that can be taken, since there are 6 unknown parameters. The simulation results are depicted in Figures 12: the estimated states and parameters converge to their true values after a transition phase due to the arbitrary choice of z(0) and w(0). The convergence speed of this observer can be tuned arbitrarily by selecting larger eigenvalues λ i 's (see [START_REF] Andrieu | Convergence speed of nonlinear luenberger observers[END_REF] for more details). From Theorem 2 and Proposition 1 it is inferred that if the input is properly selected, and if one can take a sufficiently large number of eigenvalues, the matrix P is asymptotically full rank (however, in the transient period, there is no such guarantee). This is the reason why in the algorithm, when the matrix P is too badly conditioned, the parameters are kept constant and the state estimate is obtained by integrating the model. For instance, this occurs at the beginning of the simulation (up to the time 4s). To see the observer robustness, a pseudo random noise with standard uniform distribution on the interval ±10%y(t) is added to the output signal. The figures 3 and 4 show the efficiency of the estimation algorithm with such output noise.

V. CONCLUSIONS

The design of a nonlinear Luenberger observer to estimate the state and the unknown parameters of a parametrized linear system was studied here. In a first part of the study, a Luenberger observer was shown to exist. This result is obtained from the injectivity property of a certain mapping. This result has been used on a second order linear system example and seems to offer very promising results.

Λ

  = Diag(-1, -1.25, -1.5, -1.75, -2, -2.25, -2.5) L = (1 1 1 1 1 1 1) T ; z(0) = w(0) = (0 0 0 0 0 0 0) T ( x, θ ) = ( x1 x2 θ1 θ2 θ3 θ4 )

Fig. 1 .

 1 Fig. 1. Output and its derivative: target and estimations (case without output noise).

Fig. 2 .

 2 Fig. 2. Parameters θ 1 and θ 2 : target and estimations (case without output noise).

Fig. 3 .

 3 Fig. 3. Output and its derivative: target and estimations (case with output noise).

Fig. 4 .

 4 Fig. 4. Parameters θ 3 and θ 4 : target and estimations (case with output noise).

TABLE I SYSTEM

 I

CONFIGURATION.

TABLE II OBSERVER

 II

CONFIGURATION.


V Grelet

P Dufour
email: dufour@lagep.univ-lyon1.fr

M Nadri

T Reiche

V Lemort

Modeling and control of Rankine based waste heat recovery systems for heavy duty trucks

Keywords: Waste Heat Recovery, Modeling, Control, PID. 

This paper presents a control oriented model development for waste heat recovery Rankine based control systems in heavy duty trucks. Waste heat recovery systems, such as Rankine cycle, are promising solutions to improve the fuel efficiency of heavy duty engines. Due to the highly transient operating conditions, improving the control strategy of those systems is an important step to their integration into a vehicle. The system considered here is recovering heat from both EGR and exhaust in a serial arrangement and use a mixture of water and ethanol as working fluid. The paper focuses on a comparison of a classical PID controller which is the state of the art in the automotive industry and a nonlinear model based controller in a simulation environment. The nonlinear model based controller shows better performance than the PID one and ensures safe operation of the system.

INTRODUCTION

The idea of recovering waste heat and utilize it as another form of energy is not new. Most of the actual power plants use the principle of co or even tri-generation. In the automotive industry, the most applied form of waste heat recovery system is the turbocharger which converts the waste heat into aeraulic energy by compressing the fresh air via an exhaust driven turbine. Driven by future emissions legislations and increase in fuel prices, engine gas heat recovering has recently attracted a lot of interest. In the past few years, a high number of studies have shown the interest of energy recovery Rankine based systems for heavy duty trucks engine compounding [START_REF] Sprouse | Review of Organic Rankine Cycles for Internal Combustion Engine Exhaust Waste Heat Recovery[END_REF], [START_REF] Espinosa | Contribution to the Study of Waste Heat Recovery Systems on Commercial Truck Diesel Engines[END_REF]]. Recent studies have brought a significant potential for such a system in a Heavy Duty (HD) vehicle which can lead to a decrease in fuel consumption of about 5% and reduce engine emissions. Yet many challenges have to be faced before the vehicle integration. The first challenge deals with the correct choice of fluids and system architecture [START_REF] Mago | Performance Analysis of Different Working Fluids for Use in Organic Rankine Cycles[END_REF], [START_REF] Grelet | Optimal WHR Rankine based for Heavy Duty Application[END_REF]] and shows that system simulation is a critical part of the development work. The use of water-alcohol mixture can bring some advantages in the power recuperation and overcome both disadvantages of these fluids: high freezing temperature of water and flammability of alcohol [START_REF] Latz | Comparison of Working Fluids in Both Subcritical and Supercritical Rankine Cycles for Waste Heat Recovert Systems in Heavy Duty Vehicles[END_REF]]. In those blends, Water Ethanol is quite promising and is compliant with vehicle integration where both pure fluids are not. In comparison with stationary plant, where the system is designed to run at its nominal point, the vehicle integration has to face a second challenge such as the limited cooling capacity and the highly transient behaviour of the heat sources. To deal with that, an effective control strategy is really important to maximize power recuperation and ensure a safe operation of the system. Although many papers concerning Rankine components and system optimization for mobile application can be found in the literature [START_REF] Seher | Waste Heat recovery for Commercial Vehicles with a Rankine Process[END_REF], [START_REF] Mavridou | Comparative Design Study of a Diesel Exhaust Gas Heat Exchanger for Truck Applications with Conventional and State of the Art Heat Transfer Enhancements[END_REF]], only few of them deal with control development and operating strategy improvement [START_REF] Peralez | Improving the Control Performance of an Organic Rankine Cycle System for Waste Heat Recovery from a Heavy Duty Diesel Engine using a Model Based Approach[END_REF]), Willems et al. (2012)]. One key variable to control is the working fluid temperature at the evaporator outlet / expander inlet since this temperature has a big impact on system performance [START_REF] Quoilin | Dynamic Modeling and Optimal Control Strategy of Waste Heat Recovery Organic rankine Cycles[END_REF]]. In a system perspective, it has to be as close as possible to the vapor saturation line to increase the system efficiency. An effective control of this temperature will allow to have longer recovery period by increasing the time where the expander is fed with vapour. This criterion can be achieved by reducing the standard deviation to the set point (SP). This paper is organized as follows. Section 2 presents the principle and the studied system. Section 3 approaches the modelling methodology and the resulting partial differential equation (PDE) system. Section 4 shows the two implemented controllers when section 5 compares their performances.

PRINCIPLES AND STUDIED SYSTEM

All the variables used in the following are explicitly defined in tables 2, 3 and 4 of the appendix.

Rankine Process

Rankine cycle is a widely used power generation cycle to turn heat into mechanical or electrical power. First the working fluid is pumped from a tank at the condensing pressure to the evaporator at the evaporating pressure. Then the pressurized working fluid is pre-heated, vaporized and superheated in one or several heat exchangers (HEX), also known as boilers. These HEX are linked to the heat source. The superheated vapor expands from evaporating pressure to condensing pressure in an expansion device converting the pressure and enthalpy drop into mechanical work. Finally the expanded vapour condenses through a condenser releasing heat into the heat sink (e.g. ambient air) and returns to the working fluid reservoir. In this process the changes of states in both the pump and the expander are irreversible and increase the fluid entropy to a certain extent.

Studied System

The Waste Heat Recovery System (WHRS) is compounded on a turbocharged 6 cylinder 11L 320kW HD engine using exhaust gas recirculation (EGR) and a selective catalyst reduction system (SCR) to reduce the NO x emissions. The studied Rankine cycle is recovering heat from both EGR and Exhaust applying a serial configuration of two heat exchangers. Figure 1 shows a schematic of the studied system. The mass flow rate through the two boilers is controlled by the pump speed. The expansion machine is a turbine, which has a higher power density than volumetric expanders [START_REF] Seher | Waste Heat recovery for Commercial Vehicles with a Rankine Process[END_REF], [START_REF] Lemort | A Comparison of Piston, Screw and Scroll Expanders for Small Scale Rankine Cycle Systems[END_REF]]. The working fluid is then condensed through an indirect condenser fed by coolant. Moreover the cycle is equipped with two bypass valves one located in the exhaust stream to control the amount of energy introduced in the system and a second in front of the expansion device to prevent liquid to enter in the turbine and avoid blade erosion caused by liquid droplets into a high speed rotor. The chosen working fluid is a predefined mixture of Ethanol and Water which gives the best compromise concerning performance and reduced flammability.

HEAT EXCHANGER CONTROL ORIENTED MODELING

The system dynamic is controlled by the HEX behaviour (i.e. evaporators and condenser) and models of these components are developed to dynamically predict temperature and enthalpy of transfer and working fluid at the outlet of each boiler. This is critical when coming to control design to ensure a safe operation and a proper operation of the EGR function. Safe operation means that the fluid is completely vaporized when entering the turbine in order not to destroy it. For the EGR, the aim is to have low gas temperature at the outlet of the boiler, in order not to disturb the emissions control strategy and the internal combustion process.

Model Assumptions

Several assumptions have been done to simplify the problem in a great extent. They are usually admitted when coming to heat exchanger modelling [START_REF] Feru | Experimental Validation of a Dynamic Waste Heat Recovery System Model for Control Purposes[END_REF][START_REF] Vaja | Definition of an object oriented library for the dynamic simulation of advanced energy systems: Methodologies, Tools and Applications to Combined ICE&ORC Power Plants[END_REF]]:  The transfer fluid is always considered in single phase i.e. no condensation in the EGR/exhaust gases is taken into account.  The conductive heat fluxes are neglected since the predominant phenomenon is the convection.  The pressure drops on each fluids side (transfer and working fluids) are not considered.  Both boilers are represented by a straight pipe in pipe counterflow heat exchanger, similarly to [START_REF] Vaja | Definition of an object oriented library for the dynamic simulation of advanced energy systems: Methodologies, Tools and Applications to Combined ICE&ORC Power Plants[END_REF], divided into n lumped sub-volumes in the longitudinal direction.  Fluid properties are evaluated at the outlet of each sub-volume i.e. linear profile is considered between inlet and outlet of each node.  Pressure dynamics is neglected since its time scale is very small considered to the HEX time scale.  Working fluid mass flow rate is supposed constant along the heat exchanger.

Governing Equations

Since the mass is assumed constant the continuity equation is neglected and the model is only based on the energy conservation for the working fluid (1), the gas (2) and the energy balance at the internal (3) and external wall (4).

𝜕𝑚̇𝑤 𝑓 ℎ 𝑤𝑓 𝜕𝑧 -Q ̇conv wf int = ρ 𝑤𝑓 V 𝑤𝑓 𝜕ℎ 𝑤𝑓 𝜕𝑡 .
(1)

𝜕m ̇𝑔 𝑐 𝑝 𝑔 (𝑇 𝑔 )T 𝑔 𝜕𝑧 -Q ̇conv g int -Q ̇conv g ext = ρ 𝑔 V 𝑔 𝑐 𝑝 𝑔 (𝑇 𝑔 ) ∂T 𝑔 ∂t .
(2)

Q ̇conv g int + Q ̇conv wf int = ρ wall V wall int c p wall ∂T wall int ∂t .

(3)

Q ̇conv g ext + Q ̇conv amb ext = ρ wall V wall ext c p wall ∂T wall ext ∂t , (4)
with

Q ̇conv j k = 𝛼 𝑗 𝐴 𝑒𝑥𝑐ℎ 𝑗 𝑘 (𝑇 𝑗 -𝑇 𝑤𝑎𝑙𝑙 𝑘 ),
(5) where j = g, wf, amb and k = int, ext. Furthermore, to complete the system we need boundary and initial conditions. Time-dependent boundary conditions are used in z=0 and z=L:

𝑚̇w f (𝑡, 0) = 𝑚̇w f 0 (𝑡), ℎ 𝑤𝑓 (𝑡, 0) = ℎ 𝑤𝑓 0 (𝑡), 𝑚̇g(𝑡, 𝐿) = 𝑚̇g L (𝑡), 𝑇 𝑔 (𝑡, 𝐿) = 𝑇 𝑔 𝐿 (𝑡).
The initial conditions for the gas and wall temperatures and working fluid enthalpy are given by:

𝑇 𝑔 (0, 𝑧) = 𝑇 𝑔 𝑖𝑛𝑖𝑡 (𝑧), 𝑇 𝑤𝑎𝑙𝑙 𝑖𝑛𝑡 (0, 𝑧) = 𝑇 𝑤𝑎𝑙𝑙 𝑖𝑛𝑡 𝑖𝑛𝑖𝑡 (𝑧), ℎ 𝑤𝑓 (0, 𝑧) = ℎ 𝑤𝑓 𝑖𝑛𝑖𝑡 (𝑧).
𝑇 𝑤𝑎𝑙𝑙 𝑒𝑥𝑡 (0, 𝑧) = 𝑇 𝑤𝑎𝑙𝑙 𝑒𝑥𝑡 𝑖𝑛𝑖𝑡 (𝑧),

The manipulated variable (MV) is the working fluid mass flow 𝑚̇w f 0 (𝑡), whereas the controlled variable (CV) is the working fluid enthalpy ℎ 𝑤𝑓 (𝑡, 𝐿). However the enthalpy is impossible to measure directly we have to use pressure and temperature measurement to compute it.

Heat Transfer

To model the convection from the transfer fluid to the pipe walls and from the internal pipe to the working fluid a heat transfer coefficient (α) is needed. The convection from a boundary to a moving fluid is usually represented by the dimensionless Nusselt number (Nu) which is the ratio of convective to conductive heat transfer.

𝑁𝑢 = 𝛼 * 𝑙 𝜆 , (6)
where l here represents a characteristic length and is in this case the hydraulic diameter. Numerous correlations to approach this number can be found in the literature and are usually derived from experiments, see for example [START_REF] Thome | Wolverine Tube Inc Engineering Data Book III[END_REF]. Those correlations depend on the flow regime, the number of phases and the geometry studied. In single phase, the following correlation is implemented:

𝑁𝑢 = 𝐴 𝑅𝑒 𝑛 𝑃𝑟 𝑚 , ( 7
)
where A is a constant Re and Pr are dimensionless number (respectively Reynolds and Prandtl number). By developing these two numbers (7) becomes:

𝑁𝑢 = 𝐴 ( 4𝑚π 𝑑 ℎ 𝜇 ) 𝑛 ( 𝑐 𝑝 𝜇 𝜆 ) 𝑚 . ( 8
)
Assuming the viscosity (µ) the specific heat (c p ), the heat conductivity (λ) and the characteristic length (l) constant in single phase region for the working fluid and the gas we get:

𝛼 𝑙𝑖𝑞 𝑤𝑓 = 𝛼 𝑟𝑒𝑓 𝑙𝑖𝑞 𝑤𝑓 𝑚̇𝑤 𝑓 𝑛 𝑙𝑖𝑞 𝑤𝑓 , (9a)
𝛼 𝑣𝑎𝑝 𝑤𝑓 = 𝛼 𝑟𝑒𝑓 𝑣𝑎𝑝 𝑤𝑓 𝑚̇𝑤 𝑓 𝑛 𝑣𝑎𝑝 𝑤𝑓 , (9b)
𝛼 𝑔 = 𝛼 𝑟𝑒𝑓 𝑔 𝑚̇𝑔 𝑛 𝑔 , (9c) Where the constant α ref . and the exponent n . have to be identified in liquid and vapor region for the working fluid and in single phase for the gas. In the two phase region, a more complex correlation ( 10) is used to enhance the single phase heat transfer coefficient [START_REF] Kleiber | VDI-Gesellschaft Verfahrenstechnik und Chemieingenieurwesen[END_REF]].

𝛼 2𝜑 𝑤𝑓 = 𝛼 𝑙𝑖𝑞 𝑤𝑓 {(1 -𝑞) 0,01 [(1 -𝑞) + 1,2𝑞 0,4 𝜌 𝑤𝑓 𝑙𝑖𝑞 𝜌 𝑤𝑓 𝑣𝑎𝑝 0,37 ] -2,2 + 𝑞 0,01 [ 𝛼 𝑤𝑓 𝑣𝑎𝑝 𝛼 𝑤𝑓 𝑙𝑖𝑞 (1 + 8(1 -𝑞) 0,7 𝜌 𝑤𝑓 𝑙𝑖𝑞 𝜌 𝑤𝑓 𝑣𝑎𝑝 0,67 )] -2 } -0,5 . (10)
This correlation corresponds to a tube arrangement and is practical since it creates continuity between single and two phase heat transfer coefficients and does not need transport properties such as viscosity or heat conductivity.

Working Fluid Properties

Similarly to [START_REF] Kleiber | VDI-Gesellschaft Verfahrenstechnik und Chemieingenieurwesen[END_REF], the working fluid properties are approximated using mathematical description. This allows not to rely on thermochemical database such as [START_REF] Lemmon | REFPROP Nist Standard Reference Database 23[END_REF] and creates continuity in derivative terms during the single / two phase transition. The fluid properties are only function of pressure and enthalpy.

 Temperature model:

𝑇 𝑤𝑓 = { 𝑎 𝑇𝑙𝑖𝑞 ℎ 𝑤𝑓 ² + 𝑏 𝑇𝑙𝑖𝑞 ℎ 𝑤𝑓 + 𝑐 𝑇𝑙𝑖𝑞 𝑖𝑓 ℎ 𝑤𝑓 ≤ ℎ 𝑠𝑎𝑡 𝑙𝑖𝑞 𝑇 𝑠𝑎𝑡 𝑙𝑖𝑞 + 𝑞(𝑇 𝑠𝑎𝑡 𝑣𝑎𝑝 -𝑇 𝑠𝑎𝑡 𝑙𝑖𝑞 ) 𝑖𝑓 ℎ 𝑠𝑎𝑡 𝑙𝑖𝑞 < ℎ 𝑤𝑓 < ℎ 𝑠𝑎𝑡 𝑣𝑎𝑝 𝑎 𝑇𝑣𝑎𝑝 ℎ 𝑤𝑓 ² + 𝑏 𝑇𝑣𝑎𝑝 ℎ 𝑤𝑓 + 𝑐 𝑇𝑣𝑎𝑝 𝑖𝑓 ℎ 𝑤𝑓 ≥ ℎ 𝑠𝑎𝑡 𝑣𝑎𝑝 }, (11)
where a . b . and c . are first order polynomial expressions function of pressure and q is the fluid quality defined as the quantity of vapor present in the two phase flow.

𝑞 = ℎ 𝑤𝑓 -ℎ 𝑠𝑎𝑡 𝑙𝑖𝑞 ℎ 𝑠𝑎𝑡 𝑣𝑎𝑝 -ℎ 𝑠𝑎𝑡 𝑙𝑖𝑞 . ( 12
)
The saturation temperature (T sat ) is approximated with the Wagner equation with adapted coefficient for liquid and vapor saturation [START_REF] Kleiber | VDI-Gesellschaft Verfahrenstechnik und Chemieingenieurwesen[END_REF]] and allows to make the transition between each phase.  Density model:

𝜌 𝑤𝑓 = { 𝑎 𝑑𝑙𝑖𝑞 ℎ 𝑤𝑓 ² + 𝑏 𝑑𝑙𝑖𝑞 ℎ 𝑤𝑓 + 𝑐 𝑑𝑙𝑖𝑞 𝑖𝑓 ℎ 𝑤𝑓 ≤ ℎ 𝑠𝑎𝑡 𝑙𝑖𝑞 1 𝑎 𝑑 2𝜑 ℎ 𝑤𝑓 + 𝑏 𝑑 2𝜑 𝑖𝑓 ℎ 𝑠𝑎𝑡 𝑙𝑖𝑞 < ℎ 𝑤𝑓 < ℎ 𝑠𝑎𝑡 𝑣𝑎𝑝 𝑎 𝑑𝑣𝑎𝑝 ℎ 𝑤𝑓 ² + 𝑏 𝑑𝑣𝑎𝑝 ℎ 𝑤𝑓 + 𝑐 𝑑𝑣𝑎𝑝 𝑖𝑓 ℎ 𝑤𝑓 ≥ ℎ 𝑠𝑎𝑡 𝑣𝑎𝑝 } . ( 13
)
In single phase (liq and vap) the coefficient a . b . and c . are evaluated thanks to third order polynomial function of pressure similar to:

𝑎 𝑑 = 𝑎 𝑑3 𝑃 𝑤𝑓 3 + 𝑎 𝑑2 𝑃 𝑤𝑓 2 + 𝑎 𝑑1 𝑃 𝑤𝑓 + 𝑎 𝑑0 .

(14) In two-phase region a d 2φ and 𝑏 d 2φ are approximated with:

a d 2φ = 1 a d 2φ1 𝑃 𝑤𝑓 + a d 2φ 0 , (15a)
𝑏 𝑑 2𝜑 = 1 𝑏 𝑑 2𝜑1 𝑃 𝑤𝑓 + 𝑏 𝑑 2𝜑 0 . (15b)
All coefficients are evaluated by fitting routines in Matlab using REFPROP as thermodynamic properties database [START_REF] Lemmon | REFPROP Nist Standard Reference Database 23[END_REF]].

Discretization

To implement such a model, we now have to discretize the continuous heat exchanger model with respect to space based on finite differences method. The HEX is split into n longitudinal cell where a backward Euler scheme in space is applied. The system of equations defining the response of the ith cell for transfer fluid, pipe and working fluid can be represented under the following formalism.

𝑥̇= 𝑓(𝑥, 𝑢), ( 16) where

𝑢 = [𝑚̇𝑤 𝑓 0 𝑃 𝑤𝑓 0 ℎ 𝑤𝑓 0 𝑚̇𝑔 𝐿 𝑇 𝑔 𝐿 ].
(17) The vector u contains the only MV (𝑚̇𝑤 𝑓 0 ) and four inputs disturbances: inlet gas mass flow and temperature (𝑚̇𝑔 𝐿 ) and (𝑇 𝑔 𝐿 ) and inlet working fluid pressure and enthalpy (𝑃 𝑤𝑓 0 ) and (ℎ 𝑤𝑓 0 ). And (18b)

x i = [h wf i T wall int i T g i T wall ext i], (18a)

CONTROL DEVELOPMENT

In this paper we focus on the temperature control with only regards to power generation. Our goal here is to improve the temperature stability around a given set point. For that we can manipulate the working fluid mass flow to regulate the working fluid enthalpy (or temperature) at the outlet of the exhaust boiler. Due to time issue, at the moment where this paper is written, no experimental tests of the control system are possible. The validation is based on the representative simulation model. The comparison is done on the tracking error either on a step load change or on a complete driving cycle. On each case, we start from the same equilibrium.

PI Controller

The main problem when considering PI or PID control structure is that identification for heat exchangers is very sensitive to the inputs considered [START_REF] Peralez | Improving the Control Performance of an Organic Rankine Cycle System for Waste Heat Recovery from a Heavy Duty Diesel Engine using a Model Based Approach[END_REF], [START_REF] Horst | Dynamic Heat Exchanger Model for Performance Prediction and Control System Design of Automotive Waste Heat recovery Systems[END_REF]]. This is well described when considering several operating points. Most of the results on PID tuning are reducing the process to a first order plus time delay (FOPTD) transfer function [START_REF] Skogestad | Simple analityc rules for model reduction and PID controller tuning[END_REF], Madhauranthakam et al. (2009)]. The dynamic relation between the variations of WF mass flow and the temperature is:

Δ𝑇 𝑤𝑓 𝐿 (𝑠) Δ𝑚̇𝑤 𝑓 0 (𝑠) = 𝐺 1 + 𝜏𝑠 𝑒 -𝐷𝑠 . (19)
For the identification a pseudo random binary sequence (PRBS) is the same for all operating points. As it can be seen in figure 2, this model structure is confirmed by comparing results of the reference model and an identified FOPTD around an operating point since on the worst point the agreement is around 94%. Figure 3 shows the FOPTD parameters for several inputs disturbance. For convenience they are represented versus the total heat flow rate entering:

𝑄 ̇𝑡𝑜𝑡 = 𝑐 𝑝 𝑔 (𝑇 * ) (𝑚̇𝑒 𝑔𝑟 (𝑇 𝑒𝑔𝑟 𝑖𝑛 -𝑇 𝑤𝑓 0 ) + 𝑚̇𝑒 𝑥ℎ (𝑇 𝑒𝑥ℎ 𝑖𝑛 -𝑇 𝑤𝑓 𝑠𝑎𝑡 )). (20)
where FOPTD parameters are evaluated on twenty different engine operating points, which correspond to the input disturbances, representing the overall engine map. As it can be seen below, the FOPTD parameters change a lot, according to the nonlinearity of the model developed in section 3. To improve the system control several PID tuning methods are compared to reduce both the tracking error and the control effort [START_REF] Skogestad | Simple analityc rules for model reduction and PID controller tuning[END_REF]), Madhuranthakam et al. (2008)]. Moreover, to obtain better results on a realistic driving cycle, gains are scheduled as a function of the working fluid mass flow rate entering in the system. The scheduling of each gain is realized by interpolating Kp, Ti and Td according to the working fluid mass flow feedback signal. A normal scheduling would be to select the appropriate gain in function of the operating point but schedule it as a function of the manipulated variable is here preferred to take into account the thermal inertia of the system [START_REF] Rasmussen | Gain Scheduled Control of an Air Conditioning System Using the Youla Parameterization[END_REF]]. Comparison on a step load change for several PID tuning method is shown in table 1. On two criteria: the integrated absolute error (IAE) and the total variation (TV) which are defined as:

IAE = ∫ |𝐶𝑉(𝑡) -𝑆𝑃(𝑡)|𝑑𝑡. ∞ 0 (

Nonlinear model inversion

Despite of their very good agreement, the FOPTD could not be inversed to be used as controller since they do not take account of the disturbances and therefore every correction action would only be done through a feedback on the manipulated variable. In order to improve the performance of the classical "PID" controller we use the inverse of the system model (16-18) as feedforward to compute the first part of the MV corresponding to the desired set point (Fig. 5). Details of nonlinear inversion are approached and detailed in [START_REF] Peralez | Improving the Control Performance of an Organic Rankine Cycle System for Waste Heat Recovery from a Heavy Duty Diesel Engine using a Model Based Approach[END_REF]]. To simplify the model inversion, the heat transfer coefficient for the WF is assumed constant in liquid and vapor phase.

𝛼 𝑤𝑓 = 𝛼 𝑤𝑓 𝑟𝑒𝑓 (𝑚̇𝑤 𝑓 𝑚𝑎𝑥 + 𝑚̇𝑤 𝑓 𝑚𝑖𝑛 ) 2 𝑛 , (23)
with maximum and minimum mass flow set as the mass flow corresponding to maximum and minimum WF pump speed: As the WF mass flow is considered homogenous in the entire heat exchanger it can be computed from the system (16). This two proposed control strategies are evaluated on a validated Rankine model (at the time of writing the experimental set up was not available). Gases data (exhaust and EGR) are coming from an 11L long haul truck on a highway driving cycle (Fig. 6). For the other perturbations, namely called h wf 0 and P wf 0 in ( 19), the EGR boiler inlet WF enthalpy is kept constant whereas the inlet WF pressure depends on the expansion turbine inlet conditions. Sensors and actuators dynamics are represented by first order models where time constant are fitted thanks to manufacturer data. As previously said, the controlled variable in the model is the WF enthalpy but as this one is not physically measurable and for convenience we track the WF temperature at the second boiler outlet. In order to well assess the performance, the set point changes during the simulation. (2008) and the nonlinear control presented in 4.2. Except at the beginning, where the PID is not handling the initial SP variation and the large change in disturbances appearing around 30s, it keeps the error within +/-10K, which could be acceptable when considering a volumetric expander, which are less sensitive to liquid droplets. Here a kinetic expansion machine is used and the temperature has to be as close as possible to the set point, which is better achieved by the nonlinear controller: the error remains in +/-3K. The control efforts seem similar but are slightly different (PID presents a delay compared to nonlinear controller). This paper reports control strategy development for waste heat recovery Rankine based system used on heavy duty engines. Comparison of both controllers shows a better accuracy for the inverse HEX model when comparing the tracking error. The objective to stabilize as close as possible the temperature around a given set point is achieved by using this strategy. We conclude that the nonlinear controller leads to the best performances, at the cost of developing an accurate first principle nonlinear model with all known parameters, which underline the problem of identification. Yet it still needs to be validated on test bench level to see the robustness with respect to modelling errors and parameters mismatch. However the model is compliant with classical control unit used in the automotive industry since the computational needs is low. Future work will focus on experimental validation and optimal control strategy.

ṁ𝑤 𝑓 0 = ρ wf 0 N 60 Cc 𝑝𝑢𝑚𝑝 . (24)

  Fig 1 Studied system

  Fig 2 FOPTD model identification validation

Fig

  Fig 5 New control structure proposed 5. COMPARISON OF CONTROLLERS PERFORMANCES

Fig 6

 6 Fig 6 Input disturbances Figure 7 and 8 show comparison of the two implemented control strategies: the PID based on Madhuranthakam et al.(2008) and the nonlinear control presented in 4.2. Except at the beginning, where the PID is not handling the initial SP variation and the large change in disturbances appearing around 30s, it keeps the error within +/-10K, which could be acceptable when considering a volumetric expander, which are less sensitive to liquid droplets. Here a kinetic expansion machine is used and the temperature has to be as close as possible to the set point, which is better achieved by the nonlinear controller: the error remains in +/-3K. The control efforts seem similar but are slightly different (PID presents a delay compared to nonlinear controller).

Fig

  Fig 7 Controllers performance comparison

  𝑤𝑓 𝑖-1 -ℎ 𝑤𝑓 𝑖 )-𝛼 𝑤𝑓 𝐴 𝑒𝑥𝑐ℎ 𝑤𝑓 (𝑇 𝑤𝑓 𝑖 -𝑇 𝑤𝑎𝑙𝑙 𝑖𝑛𝑡 𝑖 )) 𝜌 𝑤𝑓 𝑖 𝑉 𝑤𝑓 𝛼 𝑤𝑓 𝐴 𝑒𝑥𝑐ℎ 𝑤𝑓 (𝑇 𝑤𝑓 𝑖 -𝑇 𝑤𝑎𝑙𝑙 𝑖𝑛𝑡 𝑖 )+𝛼 𝑔 𝐴 𝑒𝑥𝑐ℎ 𝑔 𝑖𝑛𝑡 (𝑇 𝑔 𝑖 -𝑇 𝑤𝑎𝑙𝑙 𝑖𝑛𝑡 𝑖 ) 𝜌 𝑤𝑎𝑙𝑙 𝑐𝑝 𝑤𝑎𝑙𝑙 𝑉 𝑤𝑎𝑙𝑙 𝑖𝑛𝑡 ( 𝑚̇𝑔 𝑐 𝑝 𝑔 (𝑇 𝑔 )(𝑇 𝑔 𝑖-1 -𝑇 𝑔 𝑖 )-𝛼 𝑔 𝐴 𝑒𝑥𝑐ℎ 𝑔 𝑖𝑛𝑡 (𝑇 𝑔 𝑖 -𝑇 𝑤𝑎𝑙𝑙 𝑖𝑛𝑡 𝑖 ) -𝛼 𝑔 𝐴 𝑒𝑥𝑐ℎ 𝑔 𝑒𝑥𝑡 (𝑇 𝑔 𝑖 -𝑇 𝑤𝑎𝑙𝑙 𝑒𝑥𝑡 𝑖 ) ) 𝜌 𝑔 𝑉 𝑔 𝑐 𝑝 𝑔 (𝑇 𝑔 ) 𝛼 𝑎𝑚𝑏 𝐴 𝑒𝑥𝑐ℎ 𝑎𝑚𝑏 (𝑇 𝑎𝑚𝑏𝑖 -𝑇 𝑤𝑎𝑙𝑙 𝑒𝑥𝑡 𝑖 )+𝛼 𝑔 𝐴 𝑒𝑥𝑐ℎ 𝑔 𝑒𝑥𝑡 (𝑇 𝑔 𝑖 -𝑇 𝑤𝑎𝑙𝑙 𝑒𝑥𝑡 𝑖 ) 𝜌 𝑤𝑎𝑙𝑙 𝑐𝑝 𝑤𝑎𝑙𝑙 𝑉 𝑤𝑎𝑙𝑙 𝑒𝑥𝑡

	𝑓 𝑖 (𝑥, 𝑢) =
	(𝑚̇𝑤 𝑓 (ℎ
	[

] .

  21)

	𝑇𝑉 = ∫ 0 ∞	𝑑𝑀𝑉(𝑡) 𝑑𝑡	𝑑𝑡.	(22)
	The method proposed by Madhuranthakam et al. (2008) is
	minimizing the IAE with the highest TV, i.e. the control
	effort is more important but still acceptable. The IMC and
	Tyreus-Luyben tuning methods seem to be conservative and
	too slow since they present the highest IAE. Figure 4 shows
	the tracking error versus time for the five tuning methods.
	Tuning Method		IAE	TV
	Ziegler Nichols in Closed Loop 84.6325	0.0065
	Tyreus-Luyben		743.1427	0.0053
	[Skogestad (2003)]
	IMC		940.3763	0.0064
	SIMC		208.5099	0.0047
	[Skogestad (2003)]
	Madhuranthakam et al.		48.6754	0.0204
	[Madhuranthakam et al. (2008)]
	Table 1 PID tuning method comparison
	Fig 4 Tracking error comparison
	The PID tuning of Madhuranthakam et al. (2008) is now
	used.

ACKNOWLEDGMENT

This PhD thesis is collaboration between UCBL1, ULg and Volvo Trucks which is gratefully acknowledged for the funding. The French ministry of higher education and research for the financial support of the CIFRE PhD thesis 2012/549 is also acknowledged.