D Ghosh

X Bombois

J Huillery

G Scorletti

G Mercère

Optimal identification experiment design for LPV systems using the local approach

Keywords:

In this paper, we consider the problem of optimally designing the experimental conditions for LPV system identification with the local approach. Such an LPV system identification experiment is characterized by a number of local LTI identification experiments performed at constant values of the scheduling variable. The main contribution of this paper is to determine these constant values of the scheduling variable as well as the input spectra of the corresponding local LTI identification experiments in order to obtain a user-defined model accuracy with the least input energy.

Introduction

In this paper, we consider the problem of optimally designing the identification experiment leading to a model of a Linear Parameter Varying (LPV) system. A LPV system is a system whose parameters vary as a function of an external variable, the so-called scheduling variable [START_REF] Tòth | Modeling and identification of Linear Parameter Varying Systems[END_REF]. When the scheduling variable is kept constant, the LPV system reduces to a Linear Time-Invariant (LTI) system whose dynamics depend on the chosen value for the constant scheduling variable. A constant scheduling variable is often called an operating point in the literature (see e.g. [START_REF] Bachnas | A review on data-driven linear parameter-varying modeling approaches: A high-purity distillation column case study[END_REF]).

It is important to note that efficient control design methods have been developed for LPV systems [START_REF] Apkarian | A convex characterization of gain-scheduled H∞ controllers[END_REF][START_REF] Apkarian | Gain-scheduled H∞ control of missile via linear matrix inequalities[END_REF][START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF]. These control design methods deliver controllers whose parameters are, such as the system itself, a function of the scheduling variable. Methodologies have also been developed for the identification of such LPV systems. There exist two mainstream approaches for the identification of LPV systems, namely the local [START_REF] Wassink | LPV control for a wafer stage: beyond the theoretical solution[END_REF][START_REF] Liu | A local approach for the LPV identification of an actuated beam using piezoelectric actuators and sensors[END_REF][START_REF] Tòth | Modeling and identification of Linear Parameter Varying Systems[END_REF][START_REF] Steinbuch | Experimental modelling and LPV control of a motion system[END_REF] and the global approaches [START_REF] Bamieh | Identification of linear parameter varying models[END_REF][START_REF] Tòth | Modeling and identification of Linear Parameter Varying Systems[END_REF][START_REF] Bachnas | A review on data-driven linear parameter-varying modeling approaches: A high-purity distillation column case study[END_REF][START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF][START_REF] Gilson | A new data-based modelling method for identifying parsimonious nonlinear rainfall/flow models[END_REF]. In this paper, we will restrict attention to the so-called local approach for LPV system identification. In the first step of this local approach, the scheduling variable is kept constant at successive operating points and the corresponding local LTI models are identified using classical LTI identification [START_REF] Ljung | System Identification -Theory for the User, 2nd Edition[END_REF]. In a second step, these identified models are interpolated to deliver a model of the LPV system. This second step is generally done by determining an estimate of the time-invariant coefficients θ 0 of the (e.g. polyno-mial) function describing the dependence of the system parameters on the scheduling variable [START_REF] Tòth | Modeling and identification of Linear Parameter Varying Systems[END_REF][START_REF] Bachnas | A review on data-driven linear parameter-varying modeling approaches: A high-purity distillation column case study[END_REF].

In order to design such a LPV system identification procedure in an optimal way, we need to determine optimal values 1) for the number of local LTI identification experiments that will be performed, 2) for the operating points at which these local LTI identification experiments will be performed and 3) for the input spectra that will be used in each of these local LTI identification experiments. These experimental conditions will be determined in such a way that we can guarantee a userdefined model accuracy with the least input energy.

The main contribution of this paper is to propose a convex optimization problem to tackle this rather complex optimal experiment design problem. For this purpose, an intermediate contribution is to determine the minimal variance estimate for the parameter vector θ 0 describing the LPV system (see above), as well as an expression for the covariance matrix P θ of this estimate. The obtained expression for P θ is such that its inverse is a summation of the contributions of each of the local LTI identification experiments. This property allows one to transform classical accuracy constraints into constraints that are linear in the input spectra used in the local experiments. Consequently, similarly as for optimal experiment design in the LTI case [START_REF] Bombois | Least costly identification experiment for control[END_REF][START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF], one can use convex optimization to determine the optimal input spectra that have to be used at each operating point when the values of these operating points have been fixed a-priori. If this convex optimization problem is performed for a set of operating points that constitutes a fine grid of the scheduling space, the obtained solution will typically be sparse i.e. many of the input spectra will be identically zero. The optimal operating points among this fine grid of the scheduling space are consequently those corresponding to a nonzero input spectra and local LTI identification experiments will only be performed at those operating points.

Such an optimal LPV identification experiment design framework is novel. However, we had already tackled, in [START_REF] Khalate | Optimal experimental design for LPV identification using a local approach[END_REF][START_REF] Vizer | An H∞-norm-based approach for operating point selection and LPV model identification from local experiments[END_REF], the sub-problem of selecting optimal operating points for LPV system identification using the local approach. In [START_REF] Khalate | Optimal experimental design for LPV identification using a local approach[END_REF], this problem is tackled supposing that each local LTI identification experiment yields a perfect model of the local LTI system (i.e. the variance is supposed equal to zero). With this assumption, the optimal operating points are defined as those leading to the best approximation of a complex parameter dependence on the scheduling function by a polynomial function of low order. In [START_REF] Vizer | An H∞-norm-based approach for operating point selection and LPV model identification from local experiments[END_REF], an iterative methodology is developed to determine a set of operating points for which the corresponding local LTI systems represent a sufficient sample of the different local dynamics of the LPV system. However, unlike in the present paper, those operating points are not determined as the solution of an optimization problem and the variance of the identified local LTI models is also neglected.

Notations. The matrix

     X1 0 0 0 . . . 0 0 0 XN     
will be denoted diag(X 1 , ..., X N ) if the elements X i (i = 1, ..., N ) are scalar quantities while it will be denoted bdiag(X 1 , ..., X N ) if the elements X i (i = 1, ..., N ) are matrices.

2 Identification Procedure

True System

We consider the identification of a Single-Input Single-Output (SISO) LPV system i.e. a SISO system whose parameters vary with time as a function of an exogenous variable, the so-called scheduling variable p(t). We assume that p(t) is a scalar signal that we can manipulate for identification purpose and that can vary in a given scheduling interval [p min , p max ]. The considered SISO LPV system is as follows:

y(t) = y(t) + v(t) with A(z, ζ 0 (p(t))) y(t) = B(z, ζ 0 (p(t))) u(t) and D(z, ζ 0 (p(t))) v(t) = C(z, ζ 0 (p(t))) e(t) (1)
with u(t) the input signal, y(t) the measured output and e(t) a white noise signal of variance σ 2 e . In (1), A(z, ζ 0 ), B(z, ζ 0 ), C(z, ζ 0 ) and D(z, ζ 0 ) are polynomials in the delay operator z -1 and having arbitrary orders:

A(z, ζ 0 ) = 1 + a 0 1 z -1 + . . . + a 0 na z -na , B(z, ζ 0 ) = z -n k (b 0 0 + b 0 1 z -1 + . . . b 0 n b z -n b ), C(z, ζ 0 ) = 1 + c 0 1 z -1 + . . . + c 0 nc z -nc and D(z, ζ 0 ) = 1+d 0 1 z -1 +. . .+d 0 n d z -n d . The vector ζ 0 in (1)
ζ 0 i (p(t)) = np j=0 ζ 0 i,j f j (p(t)) i = 1, ..., n (2)
where n p is an arbitrary order and f j (p(t)) (j = 0, . . . , n p ) correspond to a set of unisolvent basis functions, e.g. f j (p(t)) = p j (t) (see e.g. [START_REF] Dankers | Informative data and identifiability for LPV-ARX prediction error identification[END_REF]). In (2), we suppose that n p is the same for each entry ζ 0 i (i = 1, ..., n) in order to simplify the notations (this is nevertheless not a necessity). Let us define by

θ 0 = (ζ 0 1,0 , ζ 0 1,1 , . . . , ζ 0 n,np )
T , the vector of dimension k ∆ = (n p + 1)n which contains all the time-invariant coefficients of these polynomial expansions. The relation between θ 0 and the time-varying parameter vector ζ 0 (p(t)) can thus be expressed as follows:

ζ 0 (p(t)) = T (p(t)) θ 0 (3)
for a given matrix T (p(t)) of dimension n × k which is uniquely function of f j (p(t)) (j = 0, ..., n p ). In the sequel, we will assume that the structure of the LPV system is known i.e. that n a , n b , n k , n c and n d are known in [START_REF] Apkarian | A convex characterization of gain-scheduled H∞ controllers[END_REF] and that n p is known in [START_REF] Apkarian | Gain-scheduled H∞ control of missile via linear matrix inequalities[END_REF]. Under this assumption, the parameter vector θ 0 entirely describes the LPV system.

LPV system identification using the local approach

We will now show how we can deduce an estimate θ of θ 0 using the local approach for LPV system identification and how we can evaluate the accuracy of this estimate θ. The first step of this local approach consists of performing a certain number (say M ) of so-called local LTI identification experiments (m = 1, ..., M ) [START_REF] Tòth | Modeling and identification of Linear Parameter Varying Systems[END_REF].

Let us describe one of these local LTI identification experiments. We first bring the scheduling variable to a given constant value (or operating point) p m . Then, while maintaining the scheduling variable at this operating point p(t) = p m , we apply an input sequence u(t) = u m (t) (t = 1, ..., N p ) of spectrum Φ um to the LPV system [START_REF] Apkarian | A convex characterization of gain-scheduled H∞ controllers[END_REF]. The duration of the experiment is thus denoted by N p . By measuring the corresponding output y(t) = y m (t), we obtain the data set Z m = {u m (t), y m (t) | t = 1, . . . , N p } 1 . Since p(t) is kept constant to p m during the local LTI experiment, the dynamics of the true system (1) can be represented by the following LTI sys-tem described by the time-invariant parameter vector ζ 0 (p m ) = T (p m )θ 0 :

y(t) = G(z, ζ 0 (p m )) u(t) + H(z, ζ 0 (p m )) e(t) (4)
with

G(z, ζ) = B(z, ζ)/A(z, ζ) and H(z, ζ) = C(z, ζ)/D(z, ζ).
For the sequel, it is important to observe that, due to the independent parametrization of the four polynomials, the vectors ζ 0 (p m ) and θ 0 can be split into a part that uniquely appears in the plant transfer function G and one that uniquely appears in the noise transfer function H. We will distinguish these two parts by the subscripts G and H, respectively. Consequently, (4) can be rewritten as

y(t) = G(z, ζ 0 G (p m )) u(t) + H(z, ζ 0 H (p m
)) e(t) with:

ζ 0 G (p m ) ζ 0 H (p m ) = T G (p m ) 0 0 T H (p m ) θ 0 G θ 0 H (5)
Based on the data set Z m and a full order model structure 4), we determine an estimate ζm of ζ 0 (p m ) using prediction error identification [START_REF] Ljung | System Identification -Theory for the User, 2nd Edition[END_REF]:

M = {G(z, ζ) H(z, ζ) | ζ ∈ R n } for the LTI rep- resentation (
ζm = arg min ζ 1 N p Np t=1 2 m (t, ζ) (6)
with m (t, ζ) = H -1 (z, ζ)(y m (t) -G(z, ζ)u m (t))
. This estimate has the following asymptotic distribution: ζm ∼ N (ζ 0 (p m ), P ζ,m ) where the covariance matrix P ζ,m > 0 can be estimated from the data Z m and ζm [START_REF] Ljung | System Identification -Theory for the User, 2nd Edition[END_REF]. Note furthermore that P ζ,m has the following theoretical expres- Ex(t) (E is the expectation operator) [START_REF] Ljung | System Identification -Theory for the User, 2nd Edition[END_REF]. Using this expression, we observe that the covariance matrix P ζ,m has a block-diagonal structure:

sion: P ζ,m = σ 2 e N Ē ψ m (t, ζ 0 (p m )) ψ T m (t, ζ 0 (p m )) -1 with ψ m (t, ζ) = -∂ m (t,ζ)
P ζ,m = bdiag(P ζ G ,m , P ζ H ,m
) and the block-diagonal elements of P ζ,m have the following expressions:

P -1 ζ G ,m = N p σ 2 e 1 2π π -π F u,m (e jω )F * u,m (e jω )Φ um (ω)dω (7)
P -1 ζ H ,m = N p 1 2π π -π
F e,m (e jω )F * e,m (e jω )dω [START_REF] Dankers | Informative data and identifiability for LPV-ARX prediction error identification[END_REF] with

F u,m (z) = 1 H(z,ζ 0 H (pm)) ∂G(z,ζ G ) ∂ζ G | ζ G =ζ 0 G (pm) and F e,m (z) = 1 H(z,ζ 0 H (pm)) ∂H(z,ζ H ) ∂ζ H | ζ H =ζ 0 H (pm)
. The above procedure is repeated for M different values of p m (m = 1, ..., M ) yielding M independent estimates ζm and their respective covariance matrix P ζ,m .

We will assume that the durations of these M local LTI experiments will be all equal to N p , but that the spectrum Φ um of the input signal can be chosen differently in these M experiments. Indeed, the LTI representation (4) will be different for each p m and it makes thus sense to use a different input spectrum for these different LTI systems. In the sequel, we will use the shorthand notations P M ∆ = {p 1 , . . . , p M } for the set of M operating points at which the local LTI identification experiments are performed and Φ(P M ) ∆ = {Φ u1 , . . . , Φ u M } for the set containing the input spectra used in these local LTI identification experiments.

The M local LTI identification experiments constitute the first step of the local approach for LPV system identification. The estimates ζm of ζ 0 (p m ) (m = 1, ..., M ) obtained in this first step will then be used in a second step to determine an estimate θ of the time-invariant parameter vector θ 0 using the mapping (3). In the literature, the estimate θ is generally deduced using an ordinary least-squares criterion [START_REF] Liu | A local approach for the LPV identification of an actuated beam using piezoelectric actuators and sensors[END_REF][START_REF] Steinbuch | Experimental modelling and LPV control of a motion system[END_REF][START_REF] Wassink | LPV control for a wafer stage: beyond the theoretical solution[END_REF]. This simple approach is not optimal since it does not take into account the relative accuracy of the different estimates ζm (i.e. the covariance matrices P ζ,m ). Since, in this paper, the accuracy of θ is a central aspect, we propose to instead use a generalized least-squares (or weighted least-squares) criterion (see the following theorem). Even though such an estimator is quite classical in estimation theory, it has, to the best of our knowledge, never been used before for the second step of the local LPV identification approach. Theorem 1 Consider expression (2) with unisolvent basis functions f j (p(t)) and the mapping [START_REF] Bachnas | A review on data-driven linear parameter-varying modeling approaches: A high-purity distillation column case study[END_REF]. Suppose that M local LTI identification experiments have been performed at the operating points P M ∆ = {p 1 , . . . , p M } yielding M independent estimates ζm of ζ 0 (p m ) = T (p m )θ 0 (m = 1, ..., M ). Suppose furthermore that these M independent estimates have the following distribution: ζm ∼ N (ζ 0 (p m ), P ζ,m ). Then, if and only if the set P M contains at least n p + 1 different operating points, we can deduce the following estimate θ of θ 0 :

θ = M m=1 T T (pm) P -1 ζ,m T (pm) -1 M m=1 T T (pm) P -1 ζ,m ζm (9)
and this estimate θ has the following properties:

(1) θ ∼ N (θ 0 , P θ ) with a strictly positive-definite covariance matrix P θ given by:

P θ = M m=1 T T (p m )P -1 ζ,m T (p m ) -1 , (10)
(2) the estimator ( 9) is the one leading to the smallest variance among all unbiased estimators using the given set of observations ζm (m = 1, ..., M ).

S ζ = ST θ 0 + δ (11)
with δ ∆ = Sδ having a distribution δ ∼ N (0, I nM ) (I nM denotes the identity matrix of dimension nM ). Due to the distribution of δ, the best linear unbiased estimator (BLUE) for θ 0 is the ordinary least-squares estimator based on the observation S ζ and the regressor ST [START_REF] Kay | Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory[END_REF][START_REF] Ljung | System Identification -Theory for the User, 2nd Edition[END_REF]. This estimator is given by

θ = T T S T ST -1 T T S T S ζ (12)
and we observe that this estimator is equivalent to [START_REF] Gilson | A new data-based modelling method for identifying parsimonious nonlinear rainfall/flow models[END_REF].

In [START_REF] Ljung | System Identification -Theory for the User, 2nd Edition[END_REF]Appendix 2], it is shown that, since δ is normally distributed, the estimate θ has also a normal distribution. By inserting [START_REF] Jansson | Input design via LMIs admitting frequency-wise model specifications in confidence regions[END_REF] in ( 12), we see that the expected value of θ is equal to θ 0 and that its covariance matrix P θ is equal to T T S T ST -1 which is equivalent to [START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF] in the statement of the theorem. The condition that P M contains at least n p + 1 different operating points p m is an important condition since it guarantees that θ can effectively be determined. Indeed, this condition is a necessary and sufficient condition for T to be full column rank when ( 2) is constructed with unisolvent basis functions. A full column rank T in turn ensures that T T S T ST = T T P -1 ζ T is strictly positive definite and thus invertible in [START_REF] Kay | Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory[END_REF]. This also equivalently ensures that P θ is strictly positive definite.

We still have to prove that the best linear unbiased estimator ( 9) is also the best unbiased estimator i.e. the estimator leading not only to the smallest variance among all linear unbiased estimators, but, in fact, among all unbiased estimators. For this purpose, we prove that the covariance matrix (10) of ( 12) corresponds to the Cramer-Rao bound which is the lower bound for the covariance matrix of all unbiased estimators. For this purpose, we recall that the Cramer-Rao bound is given by J -1 with J = E -∂ 2 ∂θ 2 log f (θ)| θ=θ 0 with E the expectation operator and f (θ) the probability density function of the observation vector S ζ. In our case, since δ ∼ N (0,

I nM ), f (θ) = 1 (2π) nM 2 e -0.5( ζ-T θ) T S T S( ζ-T θ) . Thus, J = E(T T P -1 ζ T ) = T T P -1
ζ T since T is fully deterministic. We thus effectively see that J -1 is given by [START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF] and thus that the second property given in the statement of the theorem holds.

Remark 1. For the properties of Theorem 1 to hold, N p must be chosen sufficiently large since these properties require ζm to be normally distributed and the latter is an asymptotic property in N p . Remark 2. It follows from the proof of Theorem 1 that, if the basis functions in (2) are not unisolvent, the result in Theorem 1 will still hold, but then under the condition that the considered set P M yields a full column rank matrix T ∆ = (T T (p 1 ), . . . , T T (p M )) T . In this case, choosing a set P M with n p + 1 different operating points will not necessarily be sufficient to guarantee this property of T .

The estimate θ in [START_REF] Gilson | A new data-based modelling method for identifying parsimonious nonlinear rainfall/flow models[END_REF] and its covariance matrix P θ have some interesting structural properties that will be used for the design of the optimal experimental conditions. The first structural property is the following one. Suppose that, based on M 1 local LTI identification experiments at operating points {p 1 , . . . , p M1 } (containing at least n p + 1 different values), we have obtained, via [START_REF] Gilson | A new data-based modelling method for identifying parsimonious nonlinear rainfall/flow models[END_REF], an estimate θ1 for θ 0 with a certain covariance matrix P θ1 that can be computed with [START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF]. Suppose also that, subsequently, we have performed additional local LTI identification experiments at {p M1+1 , . . . , p M1+M2 } (those operating points should not be necessarily different from the ones in {p 1 , . . . , p M1 }). Then, the expression [START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF] for the covariance matrix P θ of the estimator [START_REF] Gilson | A new data-based modelling method for identifying parsimonious nonlinear rainfall/flow models[END_REF] corresponding to the M = M 1 + M 2 operating points {p 1 , . . . , p M1 , p M1+1 , . . . , p M1+M2 } can be rewritten as:

P θ = P -1 θ1 + M1+M2 m=M1+1 T T (p m )P -1 ζ,m T (p m ) -1 (13)
where P ζ,m is the covariance matrix of the estimate ζm of the m th experiment (such as in Theorem 1).

The second property pertains to the structure of P θ (see [START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF]). First, observe that its inverse P -1 θ is the summation of a contribution of each local LTI experiment. Moreover, P θ has, like P ζ,m , a block-diagonal structure. Indeed, using [START_REF] Bombois | Design of least costly identification experiments : The main philosophy accompanied by illustrative examples[END_REF], P θ = bdiag(P θ G , P θ H ) with:

P -1 θ G = M m=1 T T G (p m )P -1 ζm,G T G (p m ) P -1 θ H = M m=1 T T H (p m )P -1 ζm,H T H (p m ) (14)
with P ζ G ,m and P ζ H ,m as defined in ( 7) and ( 8), respectively. The matrix P -1 θ,G is, via [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF], a linear function of Φ um (m = 1, ..., M ). This property will be important to formulate the optimal experiment design problem in a convex way. Remark 3. As already mentioned, in the above identification procedure, we have supposed that we know the structure of the LPV system. For this purpose, we not only have to know the order of G(z, ζ) and H(z, ζ) in ( 4), but we also have to correctly know the order n p in [START_REF] Apkarian | Gain-scheduled H∞ control of missile via linear matrix inequalities[END_REF]. Like for the choice of the orders of G(z, ζ) and H(z, ζ) [START_REF] Ljung | System Identification -Theory for the User, 2nd Edition[END_REF], the choice of n p can be validated aposteriori by testing the whiteness of the residual signals m (t, T (p m ) θ) (m = 1, ..., M ) and/or the whiteness of the residual vector S( ζ -T θ) (see the notation of ( 11)).

3 Optimal experiment design for LPV system identification using the local approach 3.1 Definition of the optimal experiment design problem

The local approach for LPV system identification presented in Section 2.2 involves a number of (important) experimental choices. Firstly, we have to determine the number M and the set of operating points P M = {p 1 , . . . , p M } at which local LTI identification experiments will be performed. Second, we have also to determine the set Φ(P M ) = {Φ u1 . . . Φ u M } of spectra that will be used in these local LTI experiments. The duration N p of the local experiment will be considered fixed a-priori (note nevertheless that N p will have to be chosen larger that a given threshold (see Section 3.2)).

In the sequel, we will develop a method in order to make these experimental choices in an optimal way. An objective of this optimal design will be to guarantee that the estimate θ satisfies some accuracy constraint(s). Like in [START_REF] Bombois | Design of least costly identification experiments : The main philosophy accompanied by illustrative examples[END_REF][START_REF] Potters | Experiment design for identification of structured linear systems[END_REF][START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF][START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF], the considered accuracy constraint(s) will be of the form

P -1 θ > R adm where R adm ≥ 0 is a given matrix of dimension k. If R adm is chosen as the diagonal matrix diag 1 σ 2 1 , . . . , 1 σ 2 k
, we impose that the standard deviation of each entry θ i of θ is smaller than σ i . We can also impose constraints on the accuracy of the local LTI system corresponding to an arbitrary operating point p. This local LTI system is described by (4) with a parameter vector ζ 0 (p) = T (p)θ 0 . The estimate of this parameter vector is ζp = T (p) θ and its covariance matrix is equal to T (p)P θ T T (p). We can e.g. impose an accuracy constraint on the frequency response of this local LTI model by slightly modifying the result in [START_REF] Bombois | Design of least costly identification experiments : The main philosophy accompanied by illustrative examples[END_REF], but a constraint on the accuracy of ζp can also be formulated as (T (p)P θ T T (p)) -1 > R adm,ζ where R adm,ζ > 0 is a given matrix of dimension n. This accuracy constraint can be rewritten, using Schur complements [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF], as

P -1 θ > R adm with R adm = T T (p)R adm,ζ T (p) ≥ 0.
Such a constraint can of course be imposed for different values of p. We may thus have multiple constraints P -1 θ > R j adm (j = 1, . . . , J). We will assume that all R j adm ≥ 0 have a block diagonal structure like P θ i.e. R j adm = bdiag(R j adm,G , R j adm,H ). Consequently, the accuracy constraints that θ will have to satisfy are:

P -1 θ G > R j adm,G j = 1, . . . , J (15)
P -1 θ H > R j adm,H j = 1, . . . , J (16)
with P θ G and P θ H as defined in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF].

We are now ready to formalize our optimal experiment design problem for the local LPV identification approach.

Optimal experiment design problem: determine the experimental conditions M , P M , and Φ(P M ) in such a way that the corresponding LPV identification experiment is the experiment that uses the least input energy J = N p M m=1 1 2π π -π Φ um (ω)dω while guaranteeing that the covariance matrix P θ of the estimate θ satisfies the accuracy constraints [START_REF] Liu | A local approach for the LPV identification of an actuated beam using piezoelectric actuators and sensors[END_REF] and [START_REF] Ljung | System Identification -Theory for the User, 2nd Edition[END_REF].

Like in the LTI case, the solution of the above optimal experiment design problem requires the knowledge of the true parameter vector θ 0 and the true noise variance σ 2 e since P θ G and P θ H in ( 15)-( 16) depend on θ 0 and σ 2 e (via ( 7)-( 8)). These unknown quantities have thus to be replaced by initial guesses obtained e.g. via an initial LPV identification experiment. We will discuss this further in Section 4.

It is also to be noted that, since R j adm = bdiag(R j adm,G , R j adm,H ) ≥ 0, the constraints ( 15) and ( 16) imply that P -1 θ is strictly positive-definite. The set P M solving the above optimal experiment design problem will therefore contain at least n p + 1 different operating points since the latter is a necessary (and sufficient) condition for [START_REF] Hjalmarsson | System identification of complex and structured systems[END_REF] to be strictly positive definite (see Theorem 1). 3.2 Determination of the optimal spectra for fixed operating points As we will see in the sequel, it is interesting to first consider the above optimal experiment design problem for the particular case where M and P M have been determined a-priori. The optimal spectra Φ opt (P M ) can then be determined as the solution of the following optimization problem: min Φ(P M ) J subject to [START_REF] Liu | A local approach for the LPV identification of an actuated beam using piezoelectric actuators and sensors[END_REF] and ( 16) [START_REF] Potters | Experiment design for identification of structured linear systems[END_REF] The properties of this optimization problem will be given in Theorem 2. Note beforehand that, in order to satisfy ( 16), the duration N p of each local LTI identification experiment must in any case satisfy:

N p > max j λ max   M m=1 X H (p m ) -1 R j adm,H   (18) with λ max (A) the largest eigenvalue of A and X H (p m ) ∆ = T T H (p m ) 1 2π π -π F e,m (e jω )F * e,m (e jω )dω T H (p m ) (m = 1, . . . , M ).
Theorem 2 Suppose that both M and the set P M = {p 1 , . . . , p M } are fixed. Suppose also that this set P M contains at least n p + 1 different operating points (see Theorem 1). Consider the optimization problem [START_REF] Potters | Experiment design for identification of structured linear systems[END_REF] and denote by J opt (M, P M ) the optimal cost of this optimization problem for the fixed M and P M . Then, a necessary and sufficient condition for [START_REF] Potters | Experiment design for identification of structured linear systems[END_REF] to have a solution is the condition [START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF]. When [START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF] holds, the optimization problem ( 17) is equivalent to the following convex optimiza-tion problem:

min Φ(P M ) J subject to (15) (19)
If ( 18) holds, the optimal cost J opt (M, P M ) of ( 17) is thus also the optimal cost of ( 19) and the value of this optimal cost does not depend on the chosen value for N p .

Proof. The constraints [START_REF] Ljung | System Identification -Theory for the User, 2nd Edition[END_REF] do not depend on the decision variables Φ um (m = 1, ..., M ) and have thus only to be satisfied. A necessary and sufficient condition for ( 16) to be satisfied is that (18) holds [START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF]. Consequently, when [START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF] holds, the optimization problem ( 17) is equivalent to [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF]. Note that both J and P -1 θ G are linear in the products N p Φ um (m = 1, ..., M ). Consequently, for fixed N p , the optimization problem ( 19) is convex. Moreover, if we make, in [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF], the changes of variables Φum = N p Φ um (m = 1, ..., M ), we obtain an optimization problem that is independent of N p . Consequently, the value of the optimal cost J opt (M, P M ) will not depend on the chosen value for N p .

The fact that the optimal cost J opt (M, P M ) does not depend on the choice of N p justifies why N p is not considered as a design variable in this paper.

In order to solve the convex optimization problem ( 19) in an efficient manner, we need a parametrization for the spectra Φ um (m = 1, . . . , M ). Here, we will use the following parametrization [START_REF] Jansson | Input design via LMIs admitting frequency-wise model specifications in confidence regions[END_REF] that corresponds to the spectrum of a signal that is generated by a white noise filtered by an arbitrary FIR filter of degree L (L is a user-chosen parameter):

Φ um (ω) = c m,0 + 2 L i=1 c m,i cos(iω) (20)
The positivity of (20) for all ω can be imposed by a Linear Matrix Inequality (LMI) constraint on the coefficients c m,i (i = 0, . . . , L) [START_REF] Jansson | Input design via LMIs admitting frequency-wise model specifications in confidence regions[END_REF]. These coefficients will become the decision variables of the optimization problem [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF]. With [START_REF] Steinbuch | Experimental modelling and LPV control of a motion system[END_REF], the objective function J in ( 19) is given by N p M m=1 c m,0 [START_REF] Bombois | Least costly identification experiment for control[END_REF][START_REF] Jansson | Input design via LMIs admitting frequency-wise model specifications in confidence regions[END_REF] which is a linear function of the decision variables. Moreover, P -1 ζ G ,m (see [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF]) can also be expressed as a linear function of the decision variables c m,i (i = 0, ..., L) [START_REF] Bombois | Least costly identification experiment for control[END_REF]. This property in turn makes of P -1 θ G (see ( 14)) a linear function of the decision variables c m,i (i = 0, ..., L) (m = 1, ..., M ). Consequently, with the parametrization [START_REF] Steinbuch | Experimental modelling and LPV control of a motion system[END_REF], the optimization problem ( 19) can be solved using LMI optimization techniques [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF].

Selection of the operating points

Let us now consider the original optimal experiment design problem of Section 3.1. For simplicity, we will assume that the local identification experiments can only be performed at the M grid points in a fine grid P grid = {p 1 , . . . , p M grid } of the continuous scheduling interval [p min , p max ]. In other words, we assume that the set P M can only be chosen as a subset of P grid . Under the above assumption and using the notation of Theorem 2, the cost J opt of the optimal LPV identification experiment is the minimal value of J opt (M, P M ) among all subsets P M of P grid that contains at least n p + 1 operating points (see Theorem 1). We have then the following result: J opt can be obtained by solving [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF] with M = M grid and with P M = P grid . In other words, J opt = J opt (M grid , P grid ). This result follows from the facts that the quantities which depend on the experiment in [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF] are J and P -1 θ G and that, for each LPV experiment with local LTI experiments at the points in a strict subset of P grid , we can find a (mathematically equivalent) experiment with local LTI experiments at all points in P grid that exactly lead to the same J and P -1 θ G . Indeed, let us consider an arbitrary experiment characterized by spectra Φ * um (m = 1, . . . , M ) and a set P M, * which is a strict subset of P grid (i.e. M < M grid ) and which contains at least n p + 1 different operating points. The equivalent experiment with local LTI experiments at all points in P grid and that leads to the same J and P -1 θ G is e.g. an experiment with input spectra equal to Φ * um for the p m that are also in P M, * and equal to 0 for the other M grid -M points.

Since J opt = J opt (M grid , P grid ), a solution of the optimal experiment design problem of Section 3.1 is an experiment characterized by M = M grid , P M = P grid and the spectra obtained as the solution of ( 19) with M = M grid and P M = P grid . This said, it is nevertheless important to observe that, due to the particular structure of the optimization problem [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF], its solution when M = M grid and P M = P grid will generally be sparse i.e. many of the spectra Φ um will be identically 0. Indeed, since J must be minimized, power will only be injected at those operating points that allow to obtain the maximum information on θ 0 . This phenomenon can presumably also be explained by the fact that the cost function J has a l 1 -norm structure and it is frequently observed that such cost functions, when minimized under convex constraints, generate a sparse solution (see e.g. [START_REF] Tropp | Just relax: Convex programming methods for identifying sparse signals in noise[END_REF]).

Suppose for further reference that the spectra Φ um of the solution of [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF] when M = M grid and P M = P grid are nonzero for the M operating points in the subset P M ∆ = {p 1 , . . . , p M} of P grid and denote by Φ( P M) the set of these M nonzero spectra ( M ≤ M grid ). Observe that local LTI identification experiments for which Φ um = 0 neither contribute to the cost J nor to the accuracy P -1 θ G . Consequently, the LPV identification experiment characterized by M , P M and Φ( P M) is also a solution of the optimal experiment design problem of Section 3.1 i.e. this experiment has also a cost J opt and satisfies the constraints [START_REF] Liu | A local approach for the LPV identification of an actuated beam using piezoelectric actuators and sensors[END_REF]. Since it will always be preferred to perform an LPV identification experiment with a small number of local LTI identification experiments and that performing local identifications with Φ um = 0 does not make sense, the LPV experiment character-ized by M , P M and Φ( P M) will be the experiment that will be deemed optimal and that will eventually be performed. Remark 4. To determine M , P M and Φ( P M), we thus solve the optimization problem [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF] with M = M grid and P M = P grid . This optimization problem has a solution for all values of N p . However, as shown in Theorem 2, this solution is also the solution of (17) only if [START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF] holds. We can thus follow the above procedure for the selection of M , P M and Φ( P M) with an initial value of N p .

Once P M determined, we can verify whether this initial N p satisfies the condition [START_REF] Potters | Experiment time minimisation under parameter accuracy constraints and time-domain signal amplitude bounds[END_REF]. If that is not the case, N p will have to be increased to a value N p,2 satisfying ( 18) and the spectra in Φ( P M) scaled with a factor Np Np,2 . Remark 5. When solving [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF] with M = M grid and P M = P grid as proposed above, it can be of interest to add the following constraints to the optimization problem: N p c m,0 ≤ α m (m = 1, ..., M grid ) with α m some user-chosen scalar constants. These constants α m > 0 (m = 1, ..., M grid ) can reflect the maximal energy that is allowed to apply at each operating point p m . Moreover, choosing a relatively small α m at operating points p m that are difficult to reach will favour the choice of other operating points. Note that, if we want to absolutely avoid some operating points, we can also remove those operating points from P grid . It is also to be noted that the additional constraints N p c m,0 ≤ α m will (generally) modify the obtained solution M , P M and Φ( P M). Remark 6. If, for a specific situation, the spectra of the solution of ( 19) with M = M grid and P M = P grid are nonzero at all M grid operating points (i.e. we obtain a non-sparse solution), it is important to stress that the corresponding LPV experiment with M grid local LTI experiments is nevertheless the solution of the optimal experiment design problem of Section 3.1 (under the assumption that P M must be chosen as a subset of P grid ).

4 Dependence of the solution on θ 0 and σ 2 e As in all experiment design problem, the optimal LPV identification experiment can only be determined using the procedure in the previous section if we know the true parameter vector θ 0 and the true noise covariance σ 2 e . Indeed, the covariance matrix P θ depends on θ 0 and σ 2 e (via ( 7)-( 8)). Since θ 0 and σ 2 e are of course unknown, they must consequently be replaced by some estimates in the expression for P θ .

Suppose for this purpose that we have performed a first LPV identification experiment using an arbitrary set P M1 and arbitrary spectra Φ(P M1 ). From the M 1 estimates ζm deduced from this initial LPV experiment, we can deduce an initial estimate σ2 e,1 of σ 2 e using the estimator σ2

e,1 = 1

M1 M1 m=1

1 Np Np t=1 2 m (t, ζm ) as well as an initial estimate θ1 of θ 0 via the estimator [START_REF] Gilson | A new data-based modelling method for identifying parsimonious nonlinear rainfall/flow models[END_REF]. We can therefore also deduce an estimate of the covariance matrix P θ1 of θ1 . If P -1 θ1 > R j adm ∀j, no need to proceed with optimal experiment design. However, if it is not the case, we will perform a second LPV experiment and we will design the local LTI identification experiments of this second LPV experiment in an optimal way. For this purpose, we will replace θ 0 by θ1 and σ 2 e by σ2 e,1 in the expression for P θ and we will replace the constraints

P -1 θ > R j adm by P -1 θ > R j adm -P -1
θ1 in the procedure of Section 3. Recall indeed that the covariance matrix P θ corresponding to the combination of the initial LPV experiment (i.e. the one yielding θ1 , σ2

e,1 and P θ1 ) and the additional LPV experiment is given by ( 13). Remark 7. The initial LPV experiment will also be generally used to determine optimal values for the orders n a , n b , n c , n d and n p of the LPV model structure (using e.g. the procedure of Remark 3). Consequently, it is good practice to distribute the operating points in P M1 over the entire scheduling space and to choose M 1 not too small. Indeed, M 1 should be larger than n p + 1 with the actual value of n p in (2), value that remains uncertain when designing the initial LPV experiment.

Numerical example

We will now illustrate our results by considering a true LPV system having the form (1) and described by the following polynomials:

A(z, ζ 0 ) = 1 + a 0 1 z -1 , B(z, ζ 0 ) = b 0 0 z -1 , C(z, ζ 0 ) = 1 and D(z, (p(t))) = 1 + d 0 1 z -1 .
We further assume that the variance σ 2 e of e(t) in ( 1) is equal to 0.5 and that we have a linear dependence of the parameters on the scheduling variable p(t) (n p = 1 and f j (p(t)) = p j (t)): a 0 1 (p(t)) = -0.9+0.1 p(t), b 0 0 (p(t)) = 10-1 p(t) and d 0 1 (p(t)) = -0.7+0.08 p(t). We have thus ζ 0 (p(t)) = [a 0 1 (p(t)), b 0 0 (p(t)), d 0 1 (p(t))] T and θ 0 = (-0.9, 0.1, 10, -1, -0.7, 0.08) T . The scheduling interval [p min p max ] is here given by [0 8].

In this example, we wish to determine an estimate θ of θ 0 with a guaranteed standard deviation for each entry of θ. We therefore choose R adm = diag 1 We suppose that the local LTI identification experiments can only be performed at the M grid = 17 In order to solve the optimal experiment design problem of Section 3.1, we first determine the spectra Φ um of an LPV experiment with M grid local LTI identification experiments by solving the optimization problem [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF] with M = M grid and P M = P grid . We indeed know (see Section 3.3) that this will deliver the optimal cost J opt . We first perform this optimization by computing the covariance matrices [START_REF] Boyd | Linear Matrix Inequalities in Systems and Control Theory[END_REF] using the exact values of θ 0 and of σ 2 e and we obtain J opt = J opt (M grid , P grid ) = 3511. Moreover, we observe that the obtained spectra Φ um (m = 1, ..., 17) are nonzero at only three values (i.e. p m = 0, 1.5 and 8) of the 17 operating points in P grid . The nonzero spectra at p m = 0, 1.5 and 8 are given in Figure 1. These spectra correspond to an input energy 18) for P M = {0, 1.5, 8}, this LPV experiment is the solution of the optimal experiment design problem defined in Section 3.1. This LPV experiment is thus the experiment with the least input energy that delivers an estimate θ whose covariance matrix P θ satisfies the constraint P -1 θ > R adm . It is interesting to note that the optimal LPV experiment is characterized by three local LTI identification experiments. Consequently, in this example, the optimal LPV experiment is characterized by more local LTI identification experiments than strictly required by the condition of Theorem 1. Indeed, according to Theorem 1, n p + 1 = 2 experiments at different p m are sufficient to obtain a strictly positive-definite P θ and thus to satisfy the accuracy constraints for some N p and Φ um . To further illustrate this, we have computed the costs J opt (M, P M ) that are obtained by solving [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF] for all possible subsets P M of P grid containing two operating points (i.e. M = 2). The minimal value of this cost is equal to 7904 and is obtained for P M = {0, 5.5}. We observe that this cost is, as expected, (much) higher than J opt = 3511. In Table 1, we also compare the optimal cost J opt = 3511 to the costs J opt (M, P M ) obtained for some other sets P M with M = 2. We observe that the costs J opt (M, P M ) can become very large in some cases. Until now in this example, we have supposed that we know θ 0 and σ 2 e in order to compute the expression of the covariance matrices. In practice, θ 0 and σ 2 e are of course unknown and we have instead to follow the procedure of Section 4. This procedure starts with an arbitrary initial LPV experiment. Here, we assume that this arbitrary LPV experiment is characterized by M 1 = 4 local LTI experiments with N p = 1000 and P M1 = {0, 3, 5.5, 8} (see Remark 7). The input spectra Φ um used at these four operating points are all equal to Φ um (ω) = 0.75 ∀ω. The cost J of this initial LPV experiment is thus equal to 3000 and this initial LPV experiment delivers an initial estimate θ1 of θ 0 , an initial estimate σ2 e,1 of σ 2 e and an initial covariance matrix P θ1 that does not satisfy the desired accuracy constraint. Consequently, it makes sense to optimally design the additional LPV experiment that is necessary to obtain an estimate θ satisfying the accuracy constraint. As proposed in Section 4, we follow almost the same procedure as in the case where θ 0 and σ 2 e known, but we now compute the covariance matrices using the approximations θ 0 ≈ θ1 , σ2

e,1 ≈ σ 2 e and we replace R adm by R adm -P -1 θ1 . The optimal second experiment will therefore depend on the obtained initial estimates θ1 , σ2

e,1 and P θ1 . In order to analyze this dependency, we have performed 500 Monte-Carlo simulations. For each Monte-Carlo simulation, we have performed the four initial local LTI experiments at P M1 = {0, 3, 5.5, 8} with different realizations of the white-noise input signals and different realizations of the noise signal e, we have deduced the estimates θ1 , σ2

e,1 and P θ1 and performed the optimal design of the second LPV experiment as presented in Section 4. We have observed that, for these Monte-Carlo simulations, the obtained set P M can vary, e.g. P M = {0, 1.5, 2, 8} or P M = {0, 1.5, 8}. The optimal cost J opt of this second LPV experiment will also vary for different Monte-Carlo simulations as shown in the histogram presented in Figure 3. The obtained values for this cost have a mean of 1741 while the maximal value is 2319. We observe that this cost is smaller than previously. The desired accuracy R adm is indeed here obtained by combining this additional LPV experiment with the initial one. In other words, this second LPV experiment is the one with the least input energy that guarantees that the estimate θ, which is determined by combining both this second LPV experiment and the initial one, has a covariance matrix satisfying P -1 θ > R adm . Let us also note that, for the worst of the 500 Monte-Carlo simulations, the total input energy used to obtain a sufficiently accurate estimate is equal to 3000 + 2319 = 5319. This energy is of course larger than the cost J opt = 3511 that we have obtained in the ideal case where θ 0 and σ 2 e are known. However, it is much smaller than the cost that is necessary to obtain P -1 θ > R adm using the experimental conditions of the initial experiment. Indeed, in order to satisfy P -1 θ > R adm using local LTI experiments at P M1 = {0, 3, 5.5, 8} with a white noise having the same power for all four experiments, we have to choose Φ um (ω) = 2.47 ∀ω (m = 1, ..., 4) when N p = 1000. This corresponds to a cost J = 9880.

Finally, we have verified that, despite the use of the initial estimate, the approach presented in Section 4 yields, in the vast majority of the Monte-Carlo simulations, an acceptable estimate θ of θ 0 . Recall for this purpose that P -1 θ > R adm implies that ( θ -θ 0 ) T R adm ( θθ 0 ) < ( θ -θ 0 ) T P -1 θ ( θ -θ 0 ) and recall also that ( θθ 0 ) T P -1 θ ( θ -θ 0 ) < 12.6 with probability 95%. Consequently, we can deem that θ is an acceptable estimate of θ 0 if ( θ -θ 0 ) T R adm ( θ -θ 0 ) < 12.6 and we have observed that this is indeed the case in 499 of the 500 performed Monte Carlo simulations.

Let us now return to the case where θ 0 and σ 2 e are known and recall that the solution is then characterized by an input energy of 1019 at p m = 0, an input energy of 1349 at p m = 1.5 and an input energy of 1143 at p m = 8. Let us now suppose that, for security reasons, the input energy cannot be larger than 1000 at any operating point. In this case, we have to consider the extra constraints described in Remark 6 with α m = 1000 (m = 1, ..., 17) and to solve the optimization problem [START_REF] Scorletti | Improved LMI conditions for gain scheduling and related problems[END_REF] with M = M grid and P M = P grid with these extra constraints. The obtained spectra Φ um (m = 1, ..., 17) are then nonzero at five values (i.e. p m = 0, 1.5, 2, 7.5 and 8) and the corresponding input energies at those five operating points are respectively 1000, 1000, 334, 209 and 1000. As expected, the total input energy J opt = 3543 is larger than the original one (i.e. 3511). It is also interesting to note that the maximum allowed energy is used at the operating points that had been selected in the original setting (i.e. the one without energy constraints).

Let us now also consider another situation i.e. the situation where the operating points p m ≥ 5 are much more difficult to reach. In this case, the original solution with an experiment at p m = 8 can be problematic. Let us therefore solve the optimization problem by removing from P grid the operating points p m ≥ 5. With this alternative setting, we obtain P M = {0, 1, 4.5}. We thus avoid the operating point p m = 8 that is more difficult to reach, but the price to pay is an important increase in the cost i.e. J opt = 7357. Consequently, the operating point p m = 8 seems quite important in this example to obtain the desired accuracy with low input energy.

Concluding remarks

This paper presents the first results of optimal experiment design for the identification of LPV systems. In this paper, we have restricted attention to the local approach for LPV system identification where a single scheduling variable is maintained constant at different operating points. Future work will consider the case of multiple scheduling variables. Indeed, if the approach of Section 3.3 based on a gridding of the scheduling space can in theory also be applied for that case, the complexity will increase and it is therefore interesting to investigate whether other approaches could be developed for the selection of the optimal operating points. We also wish to investigate optimal experiment design for the global approach for LPV system identification. Indeed, the particular choice for the variations of the scheduling variable(s) in the local approach (i.e. successive constant values) is not guaranteed to be the optimal one in order to obtain the most information on the time-invariant parameters θ 0 describing the LPV system. We will also continue our work towards accuracy constraints inspired by LPV control objectives.

  is a column vector of dimension n ∆ = n a +n b +1+n c +n d containing the coefficients of the four polynomials describing the LPV system i.e. a 0 1 , ..., b 0 0 ..., d 0 n d . All the entries ζ 0 i (i = 1, ..., n) of ζ 0 depend on the scheduling variable p(t) as follows:

1 ζ=

 1 T (p 1 ), . . . , T T (p M )) T , we can write that ζ ∼ N (T θ 0 , P ζ ) with P ζ a strictly positive definite block-diagonal matrix given by bdiag(P ζ,1 , . . . , P ζ,M ). We can also equivalently write that ζ = T θ 0 + δ with δ ∼ N (0, P ζ ). By multiplying this equation on both sides by the square root S of the inverse of P ζ (i.e. P -S T S), we obtain:

Fig. 1 .

 1 Fig. 1. Optimal power spectrum Φu m for pm = 0 (blue solid), for pm = 1.5 (red dashed) and for pm = 8 (black dashdot)

1

 1 = 0.003, σ 2 = 3.33 × 10 -4 , σ 3 = 0.0333, σ 4 = 0.0033, σ 5 = 0.1167 and σ 6 = 0.0133 (see Section 3.1).

Fig. 2 .Fig. 3 .

 23 Fig. 2. Modulus of the frequency response of G(z, ζ0(pm)) for pm = 0 (blue solid), for pm = 1.5 (red dashed) and for pm = 8 (black dashdot)

N p 1

 1 2π π -π Φ um (ω)dω = 1019 at p m = 0, an input energy of 1349 at p m = 1.5 and an input energy of 1143 at P M Jopt(M, P M ) P M = {0, 1.5, 8} 3511 P M = {0, 5.5} 7904 P M = {0, 8} 18388 P M = {1.5, 8} 32764 P M = {0, 1.5} 57451 Table 1 Jopt(M, P M ) for different choices of P M p m = 8. The Bode diagrams of the corresponding LTI transfer functions G(z, ζ 0 (p m )) are given in Figure 2. We observe that, since the dynamics of the LTI systems are quite different for these three values of p m , the optimal spectra Φ um are also quite different. Let us thus consider the LPV experiment with M = 3 local LTI experiments at P M = {0, 1.5, 8} and with the spectra Φ( P M) given in Figure 1. Since N p = 1000 satisfies (

The time index is set back to one at the beginning of each local LTI identification experiment.


Cyril Ray
email: cyril.ray@ecole-navale.fr

Arnaud Grancher
email: arnaud.grancher@ecole-navale.fr

Rémy Thibaud
email: remy.thibaud@ecole-navale.fr

Laurent Etienne
email: laurent.etienne@ecole-navale.fr

Spatio-Temporal Rule-based Analysis of Maritime Traffic

The maritime environment is diverse, open but partly ruled, and covers a large spectrum of ships from small sailing boats to tankers which generally exhibit type-related behaviours. Several real-time positioning systems have been designed for keeping track of ships movements. However data provided by these reporting systems are underutilized to help operators in charge of maritime safety and security. This paper introduces an information system relying on user-defined navigation rules and an inference engine to analyse maritime traffic and report abnormalities and suspicious behaviours to operators. The proposed approach improves existing ones by combining a spatio-temporal database and a rule-based inference engine.

INTRODUCTION

The maritime environment is a place where human activities (fishing, goods transportation...) continuously evolve and where traffic increases regularly. This ever increasing traffic leads to navigation difficulties and risks in coastal and crowded areas where numerous ships exhibit different movement objectives which can be conflicting (eg. sailing vs. fishing). The disasters and damages caused in the event of sea collisions can pose serious threats to the environment and human lives. The sea surveillance has therefore become a major concern in order to prevent disasters and damages but also to provide an answer to illegal immigration, drug smuggling, and terrorism. Consideration of this control issue by the International Maritime Organization (IMO) has partly evolved in the last decade from education, navigational rules to technical answers for traffic monitoring. Nowadays, ships are fitted out with real-time position report systems whose objective is to identify and locate vessels at distance (Automatic Identification System (AIS) for example). Combined with radar surveillance and others, real-time localisation of ships in coastal area (around 40 miles) is effective. A real-time maritime surveillance system relies on tools for extracting several levels of information gathered from different sources (sensors, cartographic maps, databases) as well as a set of rules [START_REF] Nilsson | Extracting rules from expert operators to support situation awareness in maritime surveillance[END_REF]. An example of such a rule is a predetermined metric which defines the suspicious behaviour according to the deviation of a vessel from the expected (or declared) route. Any new information has the potential of generating an update on dangerous or suspicious context; as a consequence, decision supports must be included in a perpetual process performing data storage, management and analysis. To describe rules, a surveillance system needs to include mechanisms from which operators can define the contextual situations he wants to be detected as suspicious, dangerous or abnormal. This mechanism can rely on a rule-based engine approach allowing to formalise rules and to ensure the link between the conceptual specification of a situation and its implementation [START_REF] Vandecasteele | Using Spatial Ontologies for Detecting Abnormal Maritime Behaviour[END_REF], [START_REF] Idiri | The automatic identification system of maritime accident risk using rule-based reasoning[END_REF], [START_REF] Morel | Integrated System for Interoperable sensors & Information sources for Common abnormal vessel behaviour detection & Collaborative identification of threat (I2C)[END_REF]. This approach entails a challenging combination of the cartographic context, position reports from ships and behavioural analysis into the inference engine for further processing. Several studies addressed spatial ontology to describe maritime traffic and identify dangerous and/or suspicious behaviours [START_REF] Roy | Anomaly detection in the maritime domain[END_REF], [START_REF] Vandecasteele | Using Spatial Ontologies for Detecting Abnormal Maritime Behaviour[END_REF]. Related to these behaviours, researches proposed solutions for anomaly detection from supervised approaches. They used set of recorded trajectories and situations, to define a panel of typical abnormal behaviours from statistical analysis [START_REF] Devogele | Mobility Data: Modelling, Management, and Understanding[END_REF], [START_REF] Brax | Anomaly detection in the surveillance domain[END_REF], [START_REF] Riveiro | Improving, maritime anomaly detection and situation awareness through interactive visualization[END_REF]. Mentioning the abnormal behaviours typology, surveillance systems with rule-based reasoning [START_REF] Popovich | Towards a close integration of Vessel Traffic Control Systems and Geographical Information systems[END_REF], [START_REF] Idiri | The automatic identification system of maritime accident risk using rule-based reasoning[END_REF], [START_REF] Roy | Rule-based expert system for maritime anomaly detection[END_REF] have been proposed. Nevertheless, although geo-based situational contexts are described in these works, integration of spatial analysis and geobased database queries are not implemented and still requires improvements. The information system presented in this paper aims to address this improvement. It presents an approach to implement decision support system for sea surveillance, combining analysis of static and dynamic geo-maritime information with a rule-based inference engine (Drools). It follows the "situation awareness" concept defined as "the perception of the elements in the environment within a volume of time and space, the comprehension of their meaning and the projection of their status in the near future" [START_REF] Endsley | Theoretical underpinnings of situation awareness: a critical review[END_REF]. Based on a spatio-temporal database model, the system includes traffic data records (AIS) and Electronic Navigational Charts (S-57). The spatio-temporal database approach allows to manage short and long term events and to extract information and knowledge about maritime traffic from spatial-temporal analysis. The remainder of this paper is organized as follows. Section 2 motivates the proposed approach and explains the architectural design of the information system. Section 3 gives implementation details and provides examples of spatio-temporal rules. Finally, Section 4 draws some conclusions and outlines future work.

A RULE-BASED REASONING APPROACH

Different strategies have been examined to address a computer-based detection of abnormal ship behaviours. A first approach relies on the automatic characterization of ship trajectories. It uses statistical analyses (clustering, aggregation, etc.) of historical data to generate, for instance, prototypes of typical paths [START_REF] Devogele | Mobility Data: Modelling, Management, and Understanding[END_REF]. Another approach relies on a semi-automatic characterization of ship trajectories [START_REF] Riveiro | Improving, maritime anomaly detection and situation awareness through interactive visualization[END_REF] which uses "supervised machine-learning techniques" to generate and refine the typical paths. Both approaches are essentially ascendant, meaning the surveillance rules are derived from data analyses.

Motivation

The design of a highly flexible surveillance system, where the business logic might have to be re-configured by a surveillance operator, requires different approaches than conventional location-based systems. Changing the business logic of an application requires a two-step process. Firstly, a domain expert, often a non-technical person, adjusts or amends the requirements. Then, a computer programmer encodes those requirements within the existing business logic. This two-steps process is pernicious. FFirstly, the logic is hidden within the code, incapacitating the expert to verify whether his requirements have been properly understood. Secondly, frequent additions or modifications of requirements may result in the production of so-called spaghetti code, noticeable by their numerous imbrications of conditional statement and overloads of design-patterns. Rule-based reasoning appears as particularly well suited for developing expert systems [START_REF] Sottara | Adding Uncertainty to a Rete-OO Inference Engine[END_REF]]. An expert system may be defined as a computer program able to reason over the knowledge of some domain expert in order to assist in the decision making process. Rule engines such as DROOLS (JBOSS), JRules BRMS (IBM), Haley BRMS (Oracle), can help to overcome such a situation by keeping the business logic isolated from the rest of the application. Production rules are pieces of knowledge represented under a condition-action format. A condition is often referred as a pattern (possibly spatial and/or temporal). Production rules are so-called because they generate behaviours. [START_REF] Browne | JBoss Drools Business Rules[END_REF] defines production rules as any fragments of knowledge that can be expressed in the format: "WHEN conditions are verified THEN perform some actions", where the WHEN part is referred to as the "left-hand side", and the THEN part as the "right-hand side". This format allows an expert to express his knowledge in a straightforward way, without computing any Java or Python code; hence, removing the need for a computer programmer to assist the expert in encoding his knowledge.

Architecture

A computing architecture has been designed for the real-time monitoring and analysis of maritime traffic with an original combination of spatio-temporal databases and expert system. The platform developed is a Java-based software with a multi-tiers architecture, and organised through a distributed data and processing model (Figure 1). The persistence layer relies on PostgreSQL database management system for data manipulation and storage. The PostGIS spatial extension is used to define and handle geographic data types and functions (e.g. distance or intersection between geometries). This layer comprises several database schemas. One has been defined for user-defined geographic information (e.g. a line ships should not cross) and one for each official nautical chart (Fig. 1, geographic data). A schema is used for positioning data (Fig. 1, real-time data) and another one has been defined for historical positioning data (Fig. 1, historical data). The temporal extent of the content can be parameterized by the user (typically 5 to 30 min). Finally editable rules defining maritime events to be detected are packaged and stored as files. The application tier allows for real-time integration of positioning information and spatio-temporal analysis and inference of ship behaviours. It comprises three different middleware. The first middleware includes parsers dedicated to the integration of positioning data coming from the AIS and to the integration of official S-57 nautical charts. The first program is a Java software, the latter is a script based on GDAL and OGR allowing for integration of vector charts. The second middleware uses Geoserver to exchange data between the databases and the end-user (presentation layer). The third middleware relies on the Drools expert system. This rule-based system is composed of three parts: a working memory filled with facts, a knowledge base filled with rules, and a rule engine matching facts with rules conditions, a process known as pattern-matching. It uses a "Data Driven" approach, meaning that rules are triggered by injection of data into the working memory. Drools continuously evaluates rules using geospatial data, the result is a set of activated rules having their conditions fully matched. The rule-based contextual annotation and associated XML-based exchange language aim to provide meta information to the presentation layer. For instance, a ship detected in a restricted area will be notified with an XML element containing the ship id, the rule(s) that triggered the notification, and the associated parameters. The presentation layer can then accommodate the visualisation of the maritime situation accordingly. The graphic user interface (GUI) displays on a single screen all the data presently visualized on many different monitors [START_REF] Rouxel | Système interactif de détection de comportements dynamiques anormaux[END_REF]. The GUI is implemented (by Mines-Telecom Institute) so as to be operated with a touch table (Diamond Touch). Thence, behavioural surveillance rules can be edited directly on the visualization interface using gestures.

EXPERT SYSTEM IMPLEMENTATION

The main objective of the expert system (based on Drools) implemented is to continuously evaluate spatial rules using geospatial data. The following section presents processes taking place in the expert system that connect geographical information dataflows and the rule-based memory management of Drools.

Implementing the system core

Figure 2 provides a representation of the main processes taking place in the expert system developed. From a user perspective (top-down approach), the surveillance operator defines rules using specific gesture on the touch table (Fig. 2-a). Then, the rules are compiled and inserted into a knowledge base (Fig. 2-b), which produces a Rete graph for Drools.

Fig. 2. Main processes involved in the system

From a data flow perspective, the AIS data received are parsed and stored in the positioning database (Fig. 2-1). The integration of spatio-temporal data within Drools memory is a two-step process. Data are firstly retrieved from the spatio-temporal database upon regular basis, and inserted into typed lists (Fig. 2-2). The mechanism has been implemented using the Java Database Connectivity (JDBC) driver. A Data Access Object (DAO) is implemented for each objects manipulated in rules (position, trajectory, coastline…). This fetching process considers two types of objects: stationary objects (S-57 data or user-defined geometries), and dynamic objects (positions and trajectories, etc.). These two types have been differentiated based on their update rates. A DAO thread regularly asks all dynamic DAO objects to check the database for the arrival of new data (every 15 seconds for AIS data). When new records are found, they are transferred into the aforementioned typed list. In a second step, a Rule Engine Runner is started and objects are pushed into the working memory of Drools (Fig. 23).

Finally, a pattern matching is performed between the facts and rules, using the Rete graph (Fig. 234). The facts injected into the working memory propagate along the graph. Ultimately, the facts that reach a leaf activate the affiliated rule (Fig. 2345). This rule is placed into the Drools agenda. At this stage, a conflict resolution strategy is adopted by the rule engine to choose which rule to execute first (Fig. 23456).

Implementing spatial rules

Two approaches have been successfully implemented and tested to execute spatial rules with our expert system. These spatial rules rely several object classes. Objects represent dynamic (e.g. ship position or trajectory) or static inputs (IHO S-57 nautical spatial objects especially restricted area, or operator-defined geometric objects: points and polylines). The first approach relies on Java Topology Suite library (JTS) which provides many different spatial operators implementing the OpenGIS Simple Features Specification. When evaluated by the inference engine, rules "call" such spatial functions and operators. The following rule (Fig. 3) illustrates that approach using JTS.

Fig. 3. Calling a spatial operation using Java Topology Suite

The rule checks whether a boat has crossed a restricted area. The left-hand side of the rule (the WHEN part) is composed of 3 conditions. Firstly, the Resare() pattern considers all the object instances of type "Resare" inside the working memory. Secondly, the AisTrajectory() pattern considers all the object instances of type AisTrajectory inside the working memory. Finally, the eval statement evaluates a boolean. Therefore the procedure called inside the parentheses must return a boolean. Here, a call is made to the "intersects" operator of JTS Topology Suite. It checks whether a polyline ($trajectory) has intersected a polygon ($resare). The right-hand side of the rule (the THEN part) logs a message each time the left-hand side is fully matched. The salience attribute represent the rule priority. The higher the salience is, the higher is the rule priority. The second approach consists in calling spatial SQL procedures which are executed remotely into PostGIS. The following rule (Figure 4) illustrates this method. The rule checks whether a point of interest (here defined by a UDG_Point, typically a ship position) is located at less than 5km from a user defined line geometry (UDG_Line). The mechanism is pretty much the same as the one explained for the previous rule. The only difference is that the condition inside the "eval" statement calls an SQL procedure. The SQL procedure uses the PostGIS function ST_Distance_Sphere to compute the linear distance in meters between a point of interest ($poi) and a polyline ($barrier) using a spherical earth.

Fig. 4. Calling a spatial SQL procedure running into PostGIS

The current objective of our implementation has been to produce of proof-of-concept for a "spatio-temporal-Drools-based" re-configurable surveillance system. We successfully tested processes connecting geographical information data-flows with Drools and fruitfully implemented spatial rules, triggering the computation of spatial operations from inside the Drools "eval" statements. However, robustness and performance of the system have not been evaluated yet. While inference engine like Drools are not initially designed for real-time processing, our approach will have to face that challenge. The number of rules, number of objects manipulated, fetching period and complexity of spatio-temporal queries are amongst parameters to evaluate.

CONCLUSION

This paper presents a maritime surveillance decision support system. This system relies on the analysis of dynamic and static geo-maritime information using reconfigurable rules. The rules are used to encapsulate knowledge about the situation awareness. This knowledge can either come from a thorough analysis of maritime data (bottom-up approach), or from the surveillance expert (top-down approach). Our system improves other existing rules-based system in the domain of maritime surveillance by defining spatial rules. A computing architecture has been designed for the real-time monitoring and analysis of maritime traffic with an original combination of spatial-temporal databases and an expert system (Drools). Spatial rules have been designed and spatial operators successfully tested. Two approaches have been implemented to run spatial rules, one based on the JTS library and the other one based on spatial SQL procedures. Following the preliminary results, further work shall address the integration of this decision support system into an operational context. As a result, the reliability and the performance of the platform could be truly evaluated.

Fig

  Fig. 1. Architecture

Acknowledgements

The research presented in this paper was supported by the European program ITEA2 and funded by French DGCIS (Reconsurve project, http://www.reconsurve.eu). The authors gratefully acknowledge this support.