Gérard Bloch

Florian Collin

Vincent Laurain

Gilles Millerioux

Marion Gilson

Éric Gnaedinger

Jean-Luc Noizette

Frederic Sarry

Mahmoud Abdelrahim

Marharyta Aleksandrova

Asma Barbata

Fengwei Chen

Laura Ciesa

Tejaswinee Darure

Thibaud Foltzer

Quentin Grandemange

Sabra Hamza

Paul Huynh

Jean Kuchly

Dawid Machala

Sarnavi Mahesh

Tatiana Manrique Espindola

Carlo Nainer

Fetra Rasoanarivo

Mathieu Rupp

Cyrine Selma

Srinivasaren- Gan Krishnan

Marylene Espinola

Patricia Espinola

Andrea Espinoza

Cezar Ferreira

Jean Martin Kyssama

Madelin Leyva

Andrey Myagkiy

Guillermo Pizarro

Silvia Robertson

Elise Vouriot

Adriana Del

Pilar Arriola

Aliz Barbosa

Vincent Bontemps

Alejandro Borroto

Jordan Cavalier

Magali De La

Cruz Barrón

Claudia De

Melo Sanchez

Mari- Ana Diaz

Isabelle Dollander

Claudia Domjahn

Carmen Gonzalez

Estibaly Hernandez

Anne-Claire Huerta

Jannes Kinscher

Mariana Limones

Karen Mancera

Damian Mar- Cual

Roberto Rama

Natalia Rebolledo

Thibaud Sauvageon

Cristina Valencia

Robert Valverde Thanks

Keywords: continuous-time model, linear time-varying model, instrumental variable, Kalman filter

The Ph.D

spent together.

When something was not on track, my friends in Chile and some other places, have been always there to give me a hand. Many thanks to you for your help and good vibes.

Finally I would like to express my profound gratitude to my family for their love and constant support. In particular, my greatest thanks go to my parents, my brother Rodrigo and my grandfather Edison. My little brother Nicolás, who rests in peace, is always in my thoughts and in my heart. I dedicate this thesis with love to him and the rest of my family.

Résumé étendu

Les personnes qui travaillent dans le domaine des Sciences et de l'ingénierie s'intéressent à la compréhension des phénomènes naturels, et à cette fin, ils construisent habituellement des modèles qui représentent le système observé. Dans le domaine de la dynamique des fluides par exemple, l'intérêt est de comprendre comment modifier la conception d'un navire afin d'améliorer sa maniabilité. Si le navire est déjà construit, le problème pour un ingénieur pourrait être de concevoir un contrôleur pour guider automatiquement le navire. Un modèle mathématique de la dynamique du navire, prenant la forme d'un ensemble d'équations, servirait alors à cette fin [START_REF] Fossen | Marine control systems[END_REF]. Une approche pour obtenir un modèle réside dans l'application des lois fondamentales de la Physique comme les lois de Newton ou les lois de Kirchhoff. La deuxième option pour la modélisation est l'identification du système, qui consiste essentiellement à construire un modèle à partir des mesures collectées de l'entrée et de la sortie du système. Cette thèse porte sur cette dernière approche.

Identification du système

La procédure d'identification des systèmes est illustrée à par la figure 1, où nous pouvons voir les trois entités principales, qui sont :

• Le jeu de données. Il s'agit d'une collection de mesures de l'entrée u et de la sortie y du système S, et est désigné par Z N = u(t k ), y(t k ) N k=1 ), où N est le nombre d'échantillons. Notez que la sortie y est contaminée par v qui peut représenter un bruit de perturbation et/ou de mesures provenant du capteur.

• La structure de modèles. C'est un ensemble désigné par M avec des structures de modèles candidates qui pourraient représenter le système [Ljung, 1999, p. 111].

Une structure de modèle est une expression mathématique, qui relie u, y et v, et est paramétrée par un vecteur θ ∈ R n θ .

• Une méthode d'identification. Compte tenu des deux premières entités, une méthode d'identification doit être choisie afin d'estimer le vecteur de paramètres θ et, par conséquent, le modèle.

Les modèles estimés à l'aide de techniques d'identification du système peuvent être classés en modèle de type boîte grise ou boîte noire [START_REF] Ljung | System Identification. Theory for the User[END_REF]. Dans le premier, les lois fondamentales de la Physique servent à construire la structure du modèle, qui comprend des paramètres inconnus. Ces paramètres sont ensuite estimés à partir d'un jeu de données et d'une méthode d'identification. Dans le cas d'un modèle de type boîte noire, la structure et ses paramètres sont directement déterminés à partir des données mesurées et la méthode d'identification. En outre, le modèle peut être à temps discret (DT) ou continu (CT). Dans ce travail, nous nous intéressons aux modèles linéaires CT qui peuvent être de type boîte grise ou boîte noire.

III

Identification hors ligne des modèles à temps continu

Dans l'estimation hors ligne, également appelée identification non récursive, ou en bloc, l'ensemble des données doit être disponible a priori. Deux approches d'identification de modèles CT linéaires invariants dans le temps (LTI) peuvent être obtenus en utilisant une estimation non récursive résumée sur le schéma de la figure 2 [START_REF] Garnier | Direct identification of continuous-time models from sampled data: issues, basic solutions and relevance. In Identification of continuous-time models from sampled data[END_REF]:

• Approche indirecte. Tout d'abord, un modèle DT est identifié puis il est converti en CT. Dans le cas LTI, cela peut être fait dans MATLAB en utilisant la fonction d2c.

• Approche directe. Il vise à estimer le modèle CT directement à partir des données échantillonnées.

Comme le montre le schéma, les paramètres CT peuvent être convertis en paramètres physiques θp , si nécessaire. Il faut remarquer que quelle que soit l'approche retenue, directe ou indirecte, il faut faire une hypothèse sur la variation des signaux entre deux instants d'échantillonnage [START_REF] Garnier | Direct identification of continuous-time models from sampled data: issues, basic solutions and relevance. In Identification of continuous-time models from sampled data[END_REF].

IV Motivés par les avantages qu'offrent les méthodes directes (voir [START_REF] Garnier | Direct identification of continuous-time models from sampled data: issues, basic solutions and relevance. In Identification of continuous-time models from sampled data[END_REF][START_REF] Garnier | The advantages of directly identifying continuous-time transfer function models in practical applications[END_REF]), nous avons décidé de travailler dans cette thèse avec ce type d'approche.

Différentes approches directes sont disponibles dans la littérature ; elles peuvent être classées comme suit [START_REF] Garnier | Continuous-time Model identification from sampled data: Implementation issues and performance evaluation[END_REF][START_REF] Garnier | Direct identification of continuous-time models from sampled data: issues, basic solutions and relevance. In Identification of continuous-time models from sampled data[END_REF]] :

• méthodes de filtrage linéaire

• méthodes intégrales

• méthodes des fonctions de modulation • d'autres approches comme les méthodes de différence finie Dans cette étude, nous considérons les méthodes de filtrage linéaire.

Identification de modèles linéaires à paramètres variables dans le temps

L'identification des modèles CT linéaires à paramètres variables (LTV) partage les aspects méthodologiques avec le cas DT LTV. Par conséquent, après quelques commentaires sur les systèmes LTV, nous discutons des approches qui ont été proposées pour l'identification DT LTV. Nous terminons la section avec l'estimation des modèles CT LTV.

Un système linéaire à paramètres variables dans le temps peut être défini comme un système dont les propriétés varient au cours du temps. Dans la vie réelle, les processus sont généralement variables dans le temps, par exemple en raison des effets de vieillissement ou des changements dans les conditions environnementales telles que la température. Parfois, la durée d'observation du système qui nous intéresse est assez courte, ce qui permet de considérer le système comme invariant dans le temps. Cependant, ce n'est pas le cas pour les systèmes que nous traitons dans cette thèse. Les paramètres d'un système LTV peuvent évoluer selon les quatre manières suivantes [Niedźwiecki, 2000, p. 60] :

• variations constantes lentes (parfois appelés dérive de paramètres)

• variations brusques peu fréquentes

• variations selon les deux modes ci-dessus (changements lents et brusques)

• tous les autres types de variations (également appelés variations de paramètres rapides)

Dans ce travail, on suppose que les paramètres varient lentement dans le temps. En outre, nous restreignons notre attention sur les systèmes linéaires (LTV).

Identification non récursive des modèles à temps discret

Les méthodes d'identification hors ligne peuvent également être utilisées pour l'identification des systèmes lentement LTV. Une approche simple et intuitive consiste à estimer des modèles LTI locaux à l'aide de la méthode des moindres carrés pondérés sur une fenêtre glissante. Deux options sont alors possibles [START_REF] Niedźwiecki | Identification of Time-varying Processes[END_REF] :

• Les moindres carrés pondérés de façon uniforme (fenêtre glissante de type rectangulaire), impliquant que tous les poids sont égaux.

V

• Les moindres carrés pondérés de façon exponentielle, où l'influence des données anciennes est atténuée au moyen d'une pondération de type exponentiel.

L'identification des systèmes LTV peut également être effectuée à l'aide de modèles déterministes pour décrire les variations de paramètres. On peut ensuite réaliser l'estimation avec l'approche de la fenêtre glissante hors ligne. Une de ces techniques est exploitée par exemple dans [START_REF] Chan | Local polynomial modeling and variable bandwidth selection for time-varying linear systems[END_REF] via une modélisation polynomiale locale. Néanmoins, nous nous intéressons à l'identification en ligne, où il est difficile de trouver une représentation appropriée des paramètres variables à l'aide de modèles déterministes. En outre, le nombre de paramètres augmente, ce qui conduit à des valeurs estimées avec une variance plus grande.

Identification récursive de modèles à temps discret

En estimation récursive, également appelée en ligne ou en temps réel, le modèle est mis à jour chaque fois que de nouvelles mesures sont disponibles. En comparaison avec l'estimation hors ligne, les algorithmes récursifs sont en général avantageux en termes de coût de calculs.

Un grand nombre de travaux de recherche a été réalisé sur l'estimation récursive des modèles DT (voir par exemple [START_REF] Ljung | Theory and Practice of Recursive Identification[END_REF][START_REF] Goodwin | Adaptive Filtering Prediction and Control[END_REF][START_REF] Ljung | Adaptation and tracking in system identification -a survey[END_REF][START_REF] Ljung | System Identification. Theory for the User[END_REF][START_REF] Haykin | Adaptive Filter Theory[END_REF]] et les références s'y rapportant). Certains algorithmes bien connus, résumés dans [START_REF] Ljung | System Identification. Theory for the User[END_REF], sont :

• Méthode récursive des moindres carrés avec facteur d'oubli (RLS-FF). Cette approche simple est dérivée des moindres carrés pondérés de manière exponentielle. L'inconvénient de cette approche est qu'elle ne convient pas si les paramètres varient à des taux différents. Dans ce cas, un facteur d'oubli de type matriciel pourrait être utilisé comme cela a été suggéré dans [START_REF] Saelid | Adaptive controller with a vector variable forgetting factor[END_REF] (voir aussi [Niedźwiecki, 2000, p. 105]).

• Filtre Kalman (KF). Cet algorithme est obtenu en supposant que les paramètres varient comme une marche aléatoire, qui est un modèle stochastique. En tant que méthode basée sur le facteur d'oubli matriciel, l'avantage de KF est qu'il convient lorsque les paramètres varient à des vitesses différentes. La marche aléatoire est le modèle stochastique le plus simple de la classe de promenade aléatoire généralisée [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF], qui a une flexibilité accrue (voir aussi [START_REF] Kitagawa | A smoothness priors-state space modeling of time series with trend and seasonality[END_REF]).

• Méthode de la variable instrumentale récursive (RIV). En tant qu'alternative non récursive, cet algorithme est algébriquement très similaire au RLS, mais avec l'avantage qu'il possède des propriétés statistiques supérieures. Différentes variantes de RIV peuvent être considérées [START_REF] Ljung | Theory and Practice of Recursive Identification[END_REF]. Pour suivre les variations de paramètres variables dans le temps, RIV peut être utilisé avec un facteur d'oubli ou le KF comme mécanisme d'adaptation [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF].

• Méthode d'erreur de prédiction récursive (RPEM). Cet algorithme est une approche basée sur le gradient. La méthode RLS et les moindres carrés (LMS) font partie de cette famille. Notez que LMS et RLS sont tout aussi bons pour le suivi des paramètres à variabilité lente [Lindbom, 1995, p. 16].

L'estimation récursive peut également être réalisée à l'aide de modèles déterministes. Comme nous l'avons déjà mentionné, dans l'identification en ligne, il est difficile de trouver une représentation appropriée des paramètres variables dans le temps à l'aide de modèles déterministes. Néanmoins, la modélisation déterministe peut offrir de bonnes performances dans le cas de variations rapides des paramètres [START_REF] Niedźwiecki | Identification of Time-varying Processes[END_REF]. Un tel problème est abordé par exemple dans [START_REF] Niedźwiecki | Fast recursive basis function estimators for time-varying processes[END_REF] en utilisant plusieurs modèles déterministes dans un schéma d'estimation parallèle. De cette façon, la difficulté de trouver un modèle déterministe est contournée. Les modèles déterministes sont également pris en compte dans le cas de certains modèles pour les systèmes non linéaires, tels que les paramètres dépendants de l'état (voir par exemple [START_REF] Young | Identification of nonlinear stochastic systems by state dependent parameter estimation[END_REF][START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF]) et les paramètres non linéaires variants (voir par exemple [START_REF] Previdi | Identification on non-linear parametrically varying models using separable least squares[END_REF]. Enfin, il est également possible de combiner des modèles déterministes et stochastiques, comme dans le cas de la modélisation de la régression harmonique dynamique [START_REF] Young | Dynamic harmonic regression[END_REF][START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF].

Les travaux de recherche pour l'estimation récursive des modèles LTV se sont poursuivis au cours des dernières décennies. Des approches basées sur la classe de modèles de promenade aléatoire généralisée sont présentées dans [START_REF] Zheng | Recursive adaptive algorithms for fast and rapidly time-varying systems[END_REF]] en considérant un algorithme basé sur RLS et dans [START_REF] Rutström | Simplified Wiener LMS Tracking with Automatic Tuning of the Step-Size[END_REF] en considérant un algorithme basé sur LMS. Un autre sujet de recherche récent est l'estimation de modèles non paramétriques dans le domaine temporel en utilisant une approche de régularisation / bayésienne (voir [START_REF] Prando | Online identification of time-varying systems: a Bayesian approach[END_REF]).

Les méthodes d'identification récursive sont nombreuses et variées. Par conséquent, choisir une méthode appropriée n'est pas une tâche facile. Les critères pour un bon choix sont par exemple le taux de convergence, la capacité de suivi, la robustesse, l'implantation informatique et les propriétés numériques (voir plus dans [Haykin, 2014, p. 23]). Dans ce travail, nous proposerons des algorithmes basés sur RLS, KF et RIV.

Identification récursive de modèles à temps continu

Une étude récente concernant l'identification des modèles DT et CT LTV est par exemple [START_REF] Lataire | Frequency Domain Measurement and Identification of Linear, Time-Varying Systems[END_REF], où une méthode fréquentielle hors ligne a été développée. Les techniques de sous-espace appliquées de manière récursive ont également été utilisées pour identifier les systèmes multivariables (MIMO) [START_REF] Bergamasco | Continuous-time predictor-based subspace identification using Laguerre filters[END_REF].

Les méthodes d'identification DT discutées ci-dessus pourraient être appliquées pour identifier un modèle de CT en utilisant l'approche indirecte. Une estimation récursive indirecte a été utilisée par exemple dans [START_REF] Gustafsson | Robust on-line estimation[END_REF].

En ce qui concerne les approches directes, une option consiste à utiliser des méthodes intégrales (voir [START_REF] Garnier | Continuous-time Model identification from sampled data: Implementation issues and performance evaluation[END_REF]), comme celles utilisées dans [START_REF] Jiang | A recursive identification method for continuous time-varying linear systems[END_REF]. La méthode CT PEM récursive présentée dans [START_REF] Östring | Recursive identification of physical parameters in a flexible robot arm[END_REF] est également une méthode intégrale, et elle a été appliquée pour estimer les paramètres physiques d'un bras de robot flexible. Une autre classe de méthodes d'identification CT directe s'appuie sur le principe des fonctions de modulation. Une implémentation récursive utilisant des fonctions de modulation est proposée dans [Co and[START_REF] Co | Batch scheme recursive parameter estimation of continuoustime systems using the modulating function method[END_REF][START_REF] Ungarala | Time-varying system identification using modulating function and spline models with application to bio-processes[END_REF]. Une troisième classe d'approches d'identification CT directe est la méthode de filtrage linéaire. Par exemple, une approche basée sur un filtre linéaire en ligne est présentée dans [START_REF] Dimogianopoulos | Adaptive control for linear slowly time-varying systems using direct least-squares estimation[END_REF], où une estimation récursive utilisant les moindres carrés pondérés de manière exponentielle est exploitée. Les algorithmes récursifs, tels que RLS, RIV et le filtre de Kalman, peuvent être couplés avec des méthodes de filtrage linéaire [Åstrom and Wittenmark, 2008, p. 59]. Cela a été proposé par [START_REF] Isermann | Identification of Dynamic Systems. An Introduction with Applications[END_REF], où la version récursive des moindres carrés associées à la méthode des filtres de variable d'état, est utilisée pour suivre non seulement les variations lentes mais également brusques de paramètres (voir aussi [START_REF] De | Recursive estimation of the continuous-time process parameters[END_REF]). Dans ce travail, nous nous concentrons sur des algorithmes récursifs basés sur des filtres linéaires, mais uniquement pour les systèmes CT dont les paramètres varient lentement dans le temps. De plus, nous développons des méthodes de variables instrumentales (IV) en raison de leurs propriétés statistiques supérieures.

List of Abbreviations

Introduction

People working in science and engineering are interested in understanding natural phenomena, and to this end they usually build models that represent the observed system.

In the field of fluid dynamics for example, the interest is to understand how to modify the design of a ship in order to improve its maneuverability. If the ship is already built, the problem for a control engineer might be to design a controller to automatically guide the vessel. A mathematical model of the ship dynamics, consisting in a set of equations, would then serve those purposes [START_REF] Fossen | Marine control systems[END_REF]. One approach to obtain a model is from physical principles like Newton's laws or Kirchhoff's laws. The second option for modeling is system identification, which basically consists in building a model based on collected measurements of the system input/output. This thesis focuses on the latter approach.

In the following, we first give some general ideas about system identification. Afterwards in Section 1.2, we present the indirect and direct approaches to identify continuoustime linear time-invariant models in an off-line fashion. Section 1.3 is dedicated to linear time-varying models. First we describe the type of linear time-varying models in which we are interested in. Putting the attention in recursive approaches, identification methods for the discrete-time and continuous-time cases are presented in Sections 1.3.2 and 1.3.3, respectively. Then, the indirect and direct identification approaches for linear time-varying systems are illustrated through a numerical example. Practical aspects and applications in recursive identification are presented in Sections 1.4 and 1.5, respectively. Finally, we present the organization of the thesis, contributions and publications.

System identification

As we already mentioned, in system identification the goal is to estimate a mathematical model of a given system (or process), based on observations of such system. The system identification procedure is shown in Figure 1.1, where we can see the three main entities, which are:

• Data set. It is a collection of measurements of the input u and output y of the system S, and it is denoted by

Z N = {u(t k ), y(t k )} N k ,
where N is the number of samples. Notice that y is contaminated by v which can represent a disturbance and/or measurement noise coming from the sensor.

• Model set. It is a set denoted by M with candidate model structures that might represent the system [Ljung, 1999, p. 111]. A model structure is a mathematical expression, which relates u, y and v, and it is parametrized by a vector θ ∈ R n θ .

• Identification method. Given the first two entities, an identification method has Figure 1.1: Basic entities in system identification: data set, model set and identification method.

to be chosen in order to estimate the parameter vector θ, and as a consequence the model.

In practice there is always a difference between the true system S and the model set M, i.e. the system does not belong to the model set (S / ∈ M). In other words, that means that there are unmodeled dynamics. Models estimated using system identification techniques can be categorized as grayor black-box models [START_REF] Ljung | System Identification. Theory for the User[END_REF]. In the former, physical principles like Newton's laws are used to build the model structure, which includes unknown parameters. Those parameters are then estimated considering a certain data set and identification method. On the other hand, in a black-box model, the structure and its parameters are solely determined by the data set and the identification method. Furthermore, the model can be in discrete-time (DT) or continuous-time (CT). In this work we are interested in CT linear models which can be of gray box or black box type.

Off-line identification of continuous-time models

In off-line estimation, also called non-recursive, batch or en-bloc identification, the whole data set must be available for estimation. In this section we illustrate with a simple example how to obtain a CT linear time-invariant model through a non-recursive identification approach.

Problem formulation

Many physical processes can be described by a CT output error (OE) model. This motivates the use of OE models in this study. As true system, let us consider then a single-input single output (SISO) CT OE representation with input u and output y, i.e.

S A o (p)x(t) = B o (p)u(t) y(t k ) = x(t k ) + e o (t k ) (1.1)
where p is the differentiation operator; x(t) is the noise-free output; and the additive term e o (t k ) is a zero-mean DT white noise sequence. The polynomials A o (p) and B o (p) are assumed to be relatively coprime and are given by

B o (p) = b o 0 p n b + b o 1 p n b -1 + . . . + b o n b (1.2) A o (p) = p na + a o 1 p na-1 + . . . + a o na (1.3)
Figure 1.2: Indirect and direct routes to identification of a CT model.

where n a ≥ n b . The argument t k in the second equation in (1.1) denotes the sampled value of the associated variable at the kth sampling instant. Measurements of the input and output with a sampling time T s form the available data set

Z N = {u(t k ); y(t k )} N k=1 .
Then, the goal is to propose a CT model structure parametrized with θ, and to estimate θ using the given data set and a certain identification method. We will assume that we know the system structure along with the polynomial degrees n a and n b , and the intersample behavior of the input.

The indirect vs direct identification methods

The two ways in which a CT model can be obtained are illustrated in the scheme of Figure 1.2; they are [START_REF] Garnier | Direct identification of continuous-time models from sampled data: issues, basic solutions and relevance. In Identification of continuous-time models from sampled data[END_REF]:

• Indirect approach. First, a DT model is identified an then it is converted to CT.

In the LTI case, this can be done in MATLAB using the routine d2c.

• Direct approach. It aims at estimating the CT model directly.

As shown in the scheme, the CT parameters can be converted to physical parameters θp , if necessary. Notice that in both the indirect and direct approaches, we need to assume a certain intersample behavior [START_REF] Garnier | Direct identification of continuous-time models from sampled data: issues, basic solutions and relevance. In Identification of continuous-time models from sampled data[END_REF]. Motivated by the advantages that offers the direct approach over the indirect approach (see [START_REF] Garnier | Direct identification of continuous-time models from sampled data: issues, basic solutions and relevance. In Identification of continuous-time models from sampled data[END_REF][START_REF] Garnier | The advantages of directly identifying continuous-time transfer function models in practical applications[END_REF]), we have decided to work in this study with the former approach. Some of the advantages of the direct approach are the following:

• It allows us to estimate directly physical quantities.

• It is suitable for non-uniformly sampled data.

• It avoids the transformation between DT and CT models which can be difficult.

• It is not sensitive to high sampling rates.

Different direct approaches are available and they can be classified as follows [START_REF] Garnier | Continuous-time Model identification from sampled data: Implementation issues and performance evaluation[END_REF][START_REF] Garnier | Direct identification of continuous-time models from sampled data: issues, basic solutions and relevance. In Identification of continuous-time models from sampled data[END_REF]]:

• linear filter methods

• integral methods

• modulating function methods

• other approaches like finite difference methods

In this study we consider linear filter methods.

Indirect identification using a DT ARX model

A simple representation for the true system S given in (1.1) is the DT autoregressive with exogenous input (ARX) model [START_REF] Ljung | System Identification. Theory for the User[END_REF]. By assuming that we know the system structure, the order of the system and the intersample behavior of the input, we can define the corresponding DT ARX representation which is given by

M : y(t k ) + α 1 y(t k-1 ) + . . . + α nα y(t k-nα ) = β 0 u(t k ) + . . . + β n β u(t k-n β ) + e(t k ) (1.4)
where e(t k ) is assumed to be a DT white noise sequence. Model (1.4) can be written as the linear regression

y(t k ) = φ T (t k )ρ + e(t k ) (1.5)
with the regressor

φ(t k ) = -y(t k-1 ) . . . -y(t k-nα ) u(t k ) . . . u(t k-n β ) T (1.6)
and the parameter vector ρ ∈ R nρ , i.e.

ρ = α 1 . . . α nα β 0 . . . β n β T (1.7)
The corresponding predictor is given by

ŷ(t k |ρ) = φ T (t k )ρ (1.8)
The parameter vector ρ can be obtained by minimizing the prediction error

ε(t k ) = y(t k ) -ŷ(t k |ρ) = y(t k ) -φ T (t k )ρ (1.9) i.e. ρ = arg min ρ∈R nρ V (1.10)
Usually, ρ is computed by choosing V as the mean-square error (MSE) defined by

V = 1 N N k=1 [y(t k ) -φ T (t k )ρ)] 2 (1.11)
This approach is the well-known prediction error method (PEM). The advantage of the DT ARX model is that (1.10) leads to the normal equations

N k=1 φ(t k )φ T (t k ) ρ = N k=1 φ(t k )y(t k ) (1.12) which have the closed-form solution ρ = N k=1 φ(t k )φ T (t k ) -1 • N k=1 φ(t k )y(t k ) (1.13)
Equation (1.13) corresponds to the least squares estimates.

Indirect identification using a DT OE model

By assuming that we know the system structure, the order of the system and the intersample behavior of the input, we can define the corresponding DT OE representation which is given by

M x(t k ) + α 1 x(t k-1 ) + . . . + α nα x(t k-nα ) = β 0 u(t k ) + . . . + β n β u(t k-n β ) y(t k ) = x(t k ) + e(t k ) (1.14)
with e(t k ) a zero-mean DT white noise sequence. The predictor for (1.14) is defined as follows [Ljung, 1999, p. 85

] ŷ(t k |ρ) = φ T (t k , ρ)ρ (1.15)
with the parameter vector ρ defined in (1.7), and the regressor

φ(t k , ρ) = -x(t k-1 , ρ) . . . x(t k-nα , ρ) u(t k ) . . . u(t k-n β ) T (1.16)
where x(t k , ρ) is obtained as the solution of the difference equation

x(t k , ρ) + α 1 x(t k-1 , ρ) + . . . + α nα x(t k-nα , ρ) = β 0 u(t k ) + . . . + β n β u(t k-n β ) (1.17)
Note that (1.15) is non-linear due to the dependence of ρ in the regressor vector. To make the difference between (1.8) and (1.15), the latter is called pseudolinear regression.

As before, the parameters can be obtained by minimizing the prediction error. However, due to the non-linear predictor (1.15), (1.10) corresponds to a non-convex optimization problem, and it cannot be guaranteed that the numerical search will yield to a global minimum. An alternative to solve the problem is to use instrumental variable techniques as for example the simplified refined instrumental variable method (SRIV) [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF]. In that approach, the non-convex optimization is replaced by iterations of convex optimization problems.

Direct identification of a CT ARX model

We assume here that the system can be represented by a CT ARX model, i.e.

M : y (na) (t k )+a 1 y (na-1) (t k )+. . .+a na y(t k ) = b 0 u (n b ) (t k )+. . .+b n b u(t k )+e(t k ) (1.18)
where x (i) (t k ) denotes the ith time derivative of the continuous-time signal x(t) at time instant t k = kT s , and e(t k ) is a zero-mean DT white noise.

The model (1.18) cannot be directly identified because the time-derivatives are usually not available. One standard approach to solve this problem is to use linear filter methods, such as the state-variable filter (SVF) (see e.g. [START_REF] Garnier | Direct identification of continuous-time models from sampled data: issues, basic solutions and relevance. In Identification of continuous-time models from sampled data[END_REF] and the prior references therein). The approach consists in filtering the model with a SVF, which is a low-pass CT filter defined by1

F (p) = 1 (p + λ svf ) na (1.19)
Then, the filtered signals y

(i) f (t k ) and u (i)
f (t k ) are obtained. They correspond to the prefiltered time-derivatives of the input and output in the bandwidth of interest, defined via the choice of λ svf . This can be written as

y (i) f (t k ) = F (p)y (i) (t k ) i = 0, . . . , n a (1.20a) u (i) f (t k ) = F (p)u (i) (t k ) i = 0, . . . , n b (1.20b)
or equivalently as

y (i) f (t k ) = p i F (p)y(t k ) i = 0, . . . , n a (1.21a) u (i) f (t k ) = p i F (p)u(t k ) i = 0, . . . , n b (1.21b)
Applying F (p) to (1.18) we obtain the following equation (which is nearly exact, except for transients arising from initial conditions),

y (na) f (t k ) + a 1 y (na-1) f (t k ) + . . . + a na y f (t k ) = b 0 u (n b ) f (t k ) + . . . + b n b u f (t k ) + v f (t k ) (1.22) with v f (t k ) = F (p)e(t k ) (1.23)
Equation (1.22) can be rewritten as a linear regression,

y (na) f (t k ) = ϕ f T (t k )θ + v f (t k ) (1.24)
where

y (na) f (t k ) = p na F (p)y(t k ) (1.25a) ϕ f T (t k ) = F (p)ϕ T (t k ) (1.25b)
with the regression vector

ϕ T (t k ) = -y (na-1) (t k ) . . . -y(t k ) u (n b ) (t k ) . . . u(t k ) (1.26)
i.e. its filtered version is given by

ϕ f T (t k ) = -y (na-1) f (t k ) . . . -y f (t k ) u (n b ) f (t k ) . . . u f (t k ) (1.27)
The parameter vector θ ∈ R n θ is defined by

θ = a 1 . . . a na b 0 . . . b n b T (1.28)
The predictor for (1.18) is then given by ŷ(na)

f (t k |θ) = ϕ f T (t k )θ (1.29)
The parameter vector (1.28) can be estimated by the least squares method, i.e. by solving the minimizing problem θ = arg min

θ∈R n θ V (1.30) with V = 1 N N k=1 [y (na) f (t k ) -ϕ f T (t k )θ] 2 (1.31)
The closed form solution to (1.30) is

θ = N k=1 ϕ f (t k )ϕ f T (t k ) -1 • N k=1 ϕ f (t k )y (na) f (t k ) (1.32)
This estimation method is called the least squares state-variable filter (LSSVF) approach.

Direct identification of a CT OE model

We assume here that the system can be represented by a CT OE model, i.e.

M x (na) (t k ) + a 1 x (na-1) (t k ) + . . . + a na x(t k ) = b 0 u (n b ) (t k ) + . . . + b n b u(t k ) y(t k ) = x(t k ) + e(t k ) (1.33)
Analogously to the DT case, the predictor y(t k |θ) corresponds to a pseudolinear regression. However, in this case it involves the solution of the differential equation

x(na) (t k , θ) + a 1 x(na-1) (t k , θ) + . . . + a na x(t k , θ) = b 0 u (n b ) (t k ) + . . . + b n b u(t k ) (1.34)
Then, the PEM estimate given by θ = arg min

θ∈R n θ 1 N N k=1 [y(t k ) -ŷ(t k |θ)] 2 (1.35)
is again a non-convex optimization problem. Alternatives to avoid the non-convex optimization are the suboptimal instrumental variable state-variable filter (IVSVF) method and the simplified refined instrumental variable method for CT models (SRIVC) [START_REF] Young | Refined instrumental variable identification of continuous-time hybrid Box-Jenkins models[END_REF], which is optimal for CT OE models. In Chapter 2, both IVSVF and SRIVC will be recalled and then adapted to linear time-varying systems.

Identification of linear time-varying models

The purpose of this section is to show how to extend the direct approach for the estimation of CT models to the case of linear time-varying (LTV) systems. The identification of CT LTV models shares methodological aspects with the DT LTV case. Therefore, after some comments about LTV systems, we discuss approaches that have been proposed in DT LTV identification.

Linear time-varying systems

A time-varying process can be defined as a system whose properties vary in time. In real-life, processes are usually time-varying, for example, due to aging effects or changes in the environmental conditions such as the temperature. Sometimes the time scale in which we are interested in, is short enough, to consider the system as time-invariant. However, this is not the case for the systems that we deal with in this thesis.

The parameters of a time-varying process can vary in the following four ways [Niedźwiecki, 2000, p. 60]:

• Slow persistent changes (sometimes called parameter drift)

• Infrequent abrupt changes

• Mixed-mode variations (slow and abrupt changes)

• All other changes (also called fast parameter variations) Note, that there is no precise distinction between fast and slow parameter variation, although several attempts to quantify the degree of non-stationarity have been done [Niedźwiecki, 2000, p. 67]. In this work, it is assumed that the parameters change slowly in time. In addition we restrain our attention to linear processes, that is, we will focus on slowly linear time-varying (LTV) systems.

LTV systems can be considered as a special class of linear parameter-varying (LPV) systems. Contrary to an LTV system, the dynamics of an LPV process depend on an independent external signal called scheduling variable, that can be measured. A typical example of an LPV process is an aircraft, where the scheduling parameters are the speed and altitude. Notice that for a given trajectory of the scheduling variable, an LTV system is obtained. Due to the similarities between these type of systems, techniques developed for LPV systems could be extended to LTV systems and vice versa. Identification of LPV models has been recently addressed in e.g. [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF][START_REF] Laurain | Direct identification of continuous-time linear parameter-varying input/output models[END_REF].

Identification of discrete-time models

Non-recursive identification methods using a sliding window

Batch estimation can also be used for the estimation of slowly LTV systems. A simple and intuitive approach is to estimate local LTI models through the weighted least squares in a sliding window. Two possible options are the following [START_REF] Niedźwiecki | Identification of Time-varying Processes[END_REF]:

• Sliding window least squares, which considers a rectangular window, meaning that all the weights are equal.

• Exponentially weighted least squares, where old data is discarded by means of exponential weights.

The identification of time-varying systems can also be tackled using deterministic models for the parameter variations, and then performing the estimation with the batch sliding window approach. Such a technique is explored for instance in [START_REF] Chan | Local polynomial modeling and variable bandwidth selection for time-varying linear systems[END_REF] by means of local polynomial modeling. Nevertheless, we are interested in online identification, where it is hard to find an appropriate representation of the timevarying parameters using deterministic models. Additionally the number of parameters is increased, leading to estimates with larger variances.

Recursive identification methods

In recursive estimation, also called on-line or real-time identification, the model is sequentially updated every time new measurements become available. In comparison with batch estimation, recursive algorithms are in general advantageous in terms of having less computational cost. In this thesis, we focus on recursive identification.

Remark 1.1 (Recursive on-line and recursive off-line identification). Note that recursive algorithms can be applied in an off-line or on-line fashion. In the former case, iterations can be used to improve the estimates (see e.g. [START_REF] Solbrand | Recursive methods for off-line identification[END_REF][START_REF] Young | Refined instrumental variable estimation: Maximum likelihood optimization of a unified Box-Jenkins model[END_REF]). Nevertheless, recursive identification is particularly suitable for on-line applications in which the model parameters are time-varying. Thus we can talk about recursive off-line or recursive on-line estimation, although in this study we concentrate on the latter. Here, whenever the word recursive is used alone, it would then mean recursive on-line.

A large number of studies have been conducted on the recursive estimation of DT models (see e.g. [START_REF] Ljung | Theory and Practice of Recursive Identification[END_REF][START_REF] Goodwin | Adaptive Filtering Prediction and Control[END_REF][START_REF] Ljung | Adaptation and tracking in system identification -a survey[END_REF][START_REF] Ljung | System Identification. Theory for the User[END_REF][START_REF] Haykin | Adaptive Filter Theory[END_REF] and the references therein). Some wellknown algorithms, summarized in [START_REF] Ljung | System Identification. Theory for the User[END_REF], are:

• Recursive least squares method with forgetting factor (RLS-FF). This simple approach is derived from the exponentially weighted least squares. The drawback of this approach is that it is not suitable if the parameters vary at different rates.

In that case, a matrix forgetting factor could be used as it has been suggested in [START_REF] Saelid | Adaptive controller with a vector variable forgetting factor[END_REF] (see also [Niedźwiecki, 2000, p. 105]).

• Kalman filter (KF). This algorithm is obtained by assuming that the parameters vary as a random walk, which is a stochastic model. As the matrix forgetting factor based method, the advantage of KF is that it is suitable when the parameters vary at different speeds. The random walk is the simplest stochastic model from the generalized random walk class [Young, 2011, p. 94], which has enhanced flexibility (see also [START_REF] Kitagawa | A smoothness priors-state space modeling of time series with trend and seasonality[END_REF]).

• Recursive instrumental variable method (RIV). As the non-recursive counterpart, this algorithm is algebraically very similar to the RLS, but with the advantage that it has superior statistical properties. Different variants of RIV can be considered [START_REF] Ljung | Theory and Practice of Recursive Identification[END_REF]]. To track time-varying parameter variations, RIV can be used together with a forgetting factor or the KF as adaptation mechanism [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF].

• Recursive prediction error method (RPEM). This algorithm is a gradient based approach. The RLS and the least mean squares (LMS) method are part of this family. Note that LMS and RLS are equally good for tracking slowly varying parameters [Lindbom, 1995, p. 16].

Recursive estimation can be also performed using deterministic models. As we already mentioned, in on-line identification, it is hard to find an appropriate representation of the time-varying parameters using deterministic models. Nevertheless, deterministic modeling can offer a good performance for the problem of fast parameter variations [START_REF] Niedźwiecki | Identification of Time-varying Processes[END_REF]. Such a problem is addressed for instance in [START_REF] Niedźwiecki | Fast recursive basis function estimators for time-varying processes[END_REF]] using more than one deterministic model in a parallel estimation scheme; in this way, the difficulty of finding a deterministic model is circumvented. Deterministic models are also considered in the case of certain models for nonlinear systems, such as state-dependent parameter (see e.g. [START_REF] Young | Identification of nonlinear stochastic systems by state dependent parameter estimation[END_REF][START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF]) and nonlinear parameter varying (see e.g. [START_REF] Previdi | Identification on non-linear parametrically varying models using separable least squares[END_REF]) models, although in this case they are used to characterize the relationship between the parameter variations and the measured system variables on which they are dependent. Finally, it is also valid to combine deterministic and stochastic models, as in the case of dynamic harmonic regression modeling [START_REF] Young | Dynamic harmonic regression[END_REF][START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF].

Research in recursive estimation of DT time-varying models has continued in the last decades, for instance, in the estimation of models with fast parameter variations. To deal with that, methods based on the generalized random walk model class have been developed. Such approaches are presented in [START_REF] Zheng | Recursive adaptive algorithms for fast and rapidly time-varying systems[END_REF] considering a RLS based algorithm and in [START_REF] Rutström | Simplified Wiener LMS Tracking with Automatic Tuning of the Step-Size[END_REF] considering a LMS based algorithm. Another recent topic of research is the estimation of non-parametric models in the time-domain, i.e. FIR models, using a regularization/Bayesian approach (see [START_REF] Prando | Online identification of time-varying systems: a Bayesian approach[END_REF] and references therein).

The amount of recursive identification methods is extensive. Therefore, choosing an appropriate method is not an easy task. Criteria to make a good choice are for instance rate of convergence, tracking ability, robustness, computational requirements and numerical properties (see more in [Haykin, 2014, p. 23]). For the reasons given above, in this work we will propose algorithms based on RLS, KF and RIV.

General recursive identification problem

Usually, the considered models can be written as a [START_REF] Ljung | System Identification. Theory for the User[END_REF] [START_REF] Ljung | System Identification. Theory for the User[END_REF]Gunnarsson, 1990, Ljung, 1999]. In this thesis, we consider algorithms that belong to (1.38) with µ(t k ) = 1. Next we review two of them, the RLS and KF methods.

Recursive least squares method

The RLS method is derived from the exponentially weighted least squares, i.e. from the following optimization problem [START_REF] Ljung | System Identification. Theory for the User[END_REF] ρ = arg min

ρ∈R nρ N ′ k=1 β(N ′ , k) y(t k ) -ŷ(t k |ρ) 2 (1.39)
where ŷ(t k |ρ) is the predictor (1.8), which corresponds to the linear regression (1.5). β(N ′ , k) is defined in terms of the forgetting factor λ(t j ),

β(N ′ , k) = N ′ j=k+1 λ(t j ) 1 ≤ k ≤ N ′ -1 1 k = N ′ (1.40)
Here we consider a constant forgetting factor λ. From (1.39) it can be shown that the recursive least squares with forgetting factor method (RLS-FF) is given by [START_REF] Ljung | System Identification. Theory for the User[END_REF]Söderström, 1983, Ljung, 1999]

ρ(t k ) = ρ(t k-1 ) + L(t k )ε(t k ) (1.41a) ε(t k ) = y(t k ) -φ T (t k )ρ(t k-1 ) (1.41b) L(t k ) = P (t k-1 )φ(t k ) λ + φ T (t k )P (t k-1 )φ(t k ) (1.41c) P (t k ) = 1 λ P (t k-1 ) -L(t k )φ T (t k )P (t k-1 ) (1.41d)
To start RLS, the initial conditions θ(t 0 ) and P (t 0 ) have to be specified. Usually, these initial conditions are set to θ(t 0 ) = 0 and P (t 0 ) = c • I nρ , where c is a constant and I is the identity matrix of order n ρ [Ljung and Söderström, 1983, p. 21].

Kalman filter method

Let us assume that the model is the linear regression (1.5), and that the parameter variations can be modeled as a random walk model; therefore the full model becomes

M ρ(t k ) = ρ(t k-1 ) + w(t k ) y(t k ) = φ T (t k )ρ(t k ) + e(t k ) (1.42)
where w(t k ) and e(t k ) are zero-mean Gaussian noises with covariance matrix Q w and variance σ 2 e , respectively. An estimate of ρ(t k ) can be obtained from a bayesian perspective, by further assuming that ρ(t k ) has a Gaussian distribution with mean ρ 0 and covariance matrix P * 0 . The KF gives then the mean ρ(t k ) and covariance matrix P * (t k ) of the posterior distribution [START_REF] Ljung | Theory and Practice of Recursive Identification[END_REF]. The KF algorithm is given by (see e.g. [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF])

Prediction step: ρ(t k |t k-1 ) = ρ(t k-1 )
(1.43a)

P * (t k |t k-1 ) = P * (t k-1 ) + Q w (1.43b)
Correction step:

ρ(t k ) = ρ(t k |t k-1 ) + L(t k )ε(t k ) (1.43c) ε(t k ) = y(t k ) -φ T (t k )ρ(t k |t k-1 ) (1.43d) L(t k ) = P * (t k |t k-1 )φ(t k ) σ 2 e + φ T (t k )P * (t k |t k-1 )φ(t k ) (1.43e) P * (t k ) = P * (t k |t k-1 ) -L(t k )φ T (t k )P * (t k |t k-1 ) (1.43f)
The normalized KF, obtained from (1.43) by considering that

P * (t k )/σ 2 e = P (t k ) (1.44)
is given by Prediction step:

ρ(t k |t k-1 ) = ρ(t k-1 ) (1.45a) P (t k |t k-1 ) = P (t k-1 ) + Q n (1.45b)
Correction step:

ρ(t k ) = ρ(t k |t k-1 ) + L(t k )ε(t k ) (1.45c) ε(t k ) = y(t k ) -φ T (t k )ρ(t k |t k-1 ) (1.45d) L(t k ) = P (t k |t k-1 )φ(t k ) 1 + φ T (t k )P (t k |t k-1 )φ(t k ) (1.45e) P (t k ) = P (t k |t k-1 ) -L(t k )φ T (t k )P (t k |t k-1 ) (1.45f)
where Q n is the normalized covariance matrix, also called noise-variance ratio matrix, and is given by andσ 2 eo . However, for other models, the KF is an ad hoc algorithm [START_REF] Ljung | System Identification. Theory for the User[END_REF]Gunnarsson, 1990, Ljung, 1999].

Q n = Q w σ 2 e (1.

Remark 1.3

The matrix Q w is usually set as

Q w = diag{σ 2 w 1 . . . σ 2 wn ρ } (1.48)
where σ 2 w i is the expected variation of the parameter ρ i (t k ). If all the parameters vary at the same speed, then Q w = σ 2 w I nρ . On the other hand, if the parameter ρ i (t k ) is known to be constant, then σ 2 w i can be set to zero or a value close to zero.

Identification of continuous-time models

A recent study concerning identification of both DT and CT LTV models is for instance [START_REF] Lataire | Frequency Domain Measurement and Identification of Linear, Time-Varying Systems[END_REF], where an off-line frequency domain method has been developed. Subspace techniques applied in a recursive fashion have been also used to identify multipleinput multiple-output (MIMO) systems [START_REF] Bergamasco | Continuous-time predictor-based subspace identification using Laguerre filters[END_REF]. However, the DT identification methods presented above could be applied to identify a CT model using the indirect approach. Indirect recursive estimation has been used for example in [START_REF] Gustafsson | Robust on-line estimation[END_REF], although this is not a common approach, because the advantages of the direct approach for LTI systems are also valid in the LTV case.

Regarding direct approaches, one option is to use integral methods (see [START_REF] Garnier | Continuous-time Model identification from sampled data: Implementation issues and performance evaluation[END_REF]), like block-pulse functions that are used in [START_REF] Jiang | A recursive identification method for continuous time-varying linear systems[END_REF]. The CT RPEM presented in [START_REF] Östring | Recursive identification of physical parameters in a flexible robot arm[END_REF] is also an integral method, and it has been applied to obtain physical parameters of a flexible robot arm. Another class of direct CT identification methods is the modulating function. A batch scheme with a recursive implementation using modulating functions is proposed in [Co and[START_REF] Co | Batch scheme recursive parameter estimation of continuoustime systems using the modulating function method[END_REF][START_REF] Ungarala | Time-varying system identification using modulating function and spline models with application to bio-processes[END_REF] for on-line estimation. A third class of direct CT identification approaches is the linear filter methods. For instance, an online linear filter based approach is presented in [START_REF] Dimogianopoulos | Adaptive control for linear slowly time-varying systems using direct least-squares estimation[END_REF], where a batch estimation using the exponentially weighted least squares is considered. Recursive algorithms, like RLS, RIV and the Kalman filter, can be coupled with linear filter methods [Åstrom and Wittenmark, 2008, p. 59]. That is done in [START_REF] Isermann | Identification of Dynamic Systems. An Introduction with Applications[END_REF], where the recursive version of LSSVF, presented in Section 1.2.2, is used to track not only slowly but also abrupt parameter variations (see also [START_REF] De | Recursive estimation of the continuous-time process parameters[END_REF]). In this study we focus on linear filter based recursive algorithms but only for CT systems whose parameters vary slowly in time.

Indirect vs direct identification

Analogously to the batch estimation case, recursive identification can be used to identify CT LTV models either in a indirect or direct way. To illustrate the two approaches, let us consider a CT LTV OE system,

S A o (p, t)x(t) = B o (p, t)u(t) y(t k ) = x(t k ) + e o (t k ) (1.49)
where e o (t k ) is a zero-mean DT white noise sequence, and

B o (p, t) = b o 0 (t)p n b + b o 1 (t)p n b -1 + . . . + b o n b (t) (1.50) A o (p, t) = p na + a o 1 (t)p na-1 + . . . + a o na (t) (1.51)
In this Chapter, we assume that system S can be represented by a CT LTV ARX model set2 , i.e.

M : y (na) (t k ) + a 1 (t k )y (na-1) (t k ) + . . . + a na (t k )y(t k ) = b 0 (t k )u (n b ) (t k ) + . . . + b n b (t k )u(t k ) + e(t k ) (1.52)
In the indirect approach, after building the predictor, the estimates are compute with (1.38). Since the system is LTV, the conversion from DT to CT has to be done in every recursion.

In the direct approach using linear filter methods, we need to write (1.52) as a linear regression. In the LTV case, this is not obvious since the differentiation operator p is not commutative with the time-dependent variables. Therefore, a filter cannot be applied in order to obtain prefiltered time-derivatives. However, if we assume that the parameters vary slowly, we can neglect the non-commutativity problem. This issue will be formalized and discussed further in Chapter 2. In the meantime, let us state that this is possible. Analogously to the LTI case (see (1.24)), applying a filter F (p) to (1.52) we obtain the filtered time-varying linear regression

y (na) f (t k ) = ϕ f T (t k )θ(t k ) + v f (t k ) (1.53)
with y

(na) f (t k ) and ϕ f (t k ) defined in (1.25), and v f (t k ) = F (p)e(t k ). The parameter vector is now

θ(t k ) = a 1 (t k ) . . . a na (t k ) b 0 (t k ) . . . b n b (t k ) T (1.54)
Then, the corresponding prediction error is

ε(t k ) = y (na) f (t k ) -ϕ T f (t k ) θ(t k-1 ) (1.55)
Once the predictor is obtained, the CT parameters are obtained with a recursive algorithm of the general class (1.38), i.e. in this case

θ(t k ) = θ(t k-1 ) + µ(t k )L(t k )ε(t k ) (1.56)
A scheme of the indirect and direct recursive identification approaches, applied to a linear regression model, is presented in Figure 1.3. For a pseudolinear regression the procedure is quite similar, but more involved due to the computation of the predictor. Note that the prefiltering is an inherent step of the linear filter methods. In DT recursive estimation, data filtering is not compulsory, but it is sometimes recommended to cope with disturbances [Åstrom and Wittenmark, 2008, p. 466] or to enhance a certain frequency band. As the p operator, the shift operator q is also not commutative with time-dependent variables. Thus, in DT recursive estimation, if the data is filtered, the non-commutativity problem mentioned above does also exist.

Recursive least squares state-variable filter method

The direct LSSVF method presented in Section 1.2.2 can be implemented in a recursive fashion based on the methodology used in the DT case. To cope with the parameter variation we can consider the forgetting factor. Then the estimates are obtained with the following algorithm [START_REF] Isermann | Identification of Dynamic Systems. An Introduction with Applications[END_REF]:

θ(t k ) = θ(t k-1 ) + L(t k )ε(t k ) (1.57a) ε(t k ) = y (na) f (t k ) -ϕ f T (t k ) θ(t k-1 ) (1.57b) L(t k ) = P (t k-1 )ϕ f (t k ) λ + ϕ f T (t k )P (t k-1 )ϕ f (t k )
(1.57c)

P (t k ) = 1 λ P (t k-1 ) -L(t k )ϕ f T (t k )P (t k-1 ) (1.57d)
where y

(na) f (t k ) and ϕ f (t k ) are defined in (1.25). Algorithm (1.57) will be called the recursive least squares state-variable filter with forgetting factor (RLSSVF-FF) method.

The second option is to couple LSSVF with KF as an adaptation mechanism for the time-varying parameters. Then, the parameters are modeled through a random walk, i.e.

θ(t k ) = θ(t k-1 ) + w(t k ) (1.58)
The full model is given by (1.52) and (1.58), with w(t k ) and e(t k ) independent zero-mean DT Gaussian noise processes with covariance matrix Q w variance σ 2 e , respectively. Then, the algorithm, based on the DT methodology (see e.g. [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF]), is defined by Prediction step:

θ(t k |t k-1 ) = θ(t k-1 )
(1.59a)

P (t k |t k-1 ) = P (t k-1 ) + Q n (1.59b)
Correction step:

θ(t k ) = θ(t k |t k-1 ) + L(t k )ε(t k ) (1.59c) ε(t k ) = y (na) f (t k ) -ϕ T f (t k ) θ(t k |t k-1 ) (1.59d) L(t k ) = P (t k |t k-1 )ϕ f (t k ) 1 + ϕ T f (t k )P (t k |t k-1 )ϕ f (t k ) (1.59e) P (t k ) = P (t k |t k-1 ) -L(t k )ϕ T f (t k )P (t k |t k-1 ) (1.59f)
where y

(na) f (t k ) and ϕ f (t k ) are defined in (1.25).
Algorithm (1.59) will be called the recursive least squares state-variable filter based Kalman filter (RLSSVF-KF) method.

Two algorithms have been presented, RLSSVF-FF and RLSSVF-KF. In the thesis, we mainly use RLSSVF-KF; therefore, to simplify the notation we denote it just as RLSSVF.

To start RLSSVF, the initial conditions θ(t 0 ) and P (t 0 ) have to be specified. In this work, unless otherwise stated we choose θ(t 0 ) = 0 and P (t 0 ) = 10 4 • I n θ .

Remark 1.4

The cut-off frequency λ svf is a user parameter that should be chosen somewhat larger than the system bandwidth. In the LTV case and especially for systems with large variations in the bandwidth, the specification of λ svf can be critical since the system bandwidth is time-varying. Some other aspects about the SVF and the choice of λ svf are discussed in Section 3.1. Notice that in CT recursive estimation, an additional computation is needed to obtain the simulated model output, which we can use to define the output error as follows

ε y (t k ) = y(t k ) -ŷ( θ(t k ))
(1.60)

Numerical example -Similar parameter variations

We illustrate the indirect and direct approach with an example where the parameters vary at a similar rate. Let us consider the following CT LTV OE system

S p 2 + a o 1 (t)p + a o 2 (t) x(t) = b o 0 (t)u(t) y(t k ) = x(t k ) + e(t k ) (1.61)
where

a o 1 (t) = 2ζ d ω n (t) a o 2 (t) = b o 0 (t) = ω 2 n (t) (1.62)
with constant damping ζ d = 0.6 and

ω n (t) = 2(1 + 0.3 cos(2πt/4500)) (1.63)
For this second order underdamped system, the measurements have a sampling time T s = 0.1 s and the total simulation time is 9000 s. Let us define the signal-to-noise ratio (SNR) as follows SNR = 10 log P x P e (1.64)

where P x and P e denote the average power of the noise-free output and the additive noise, respectively. In this example, the DT measurement noise is a zero-mean, Gaussian noise with constant variance 0.05, that corresponds to SNR ≈ 11 dB. The input of the system is a PRBS. Part of the input-output data is shown in Figure 1.4. The step responses and the Bode plots for some frozen systems obtained every 250 s are plotted in Figures 1.5 and 1.6, respectively. The bandwidth of the same frozen system are shown in Figure 1.7.

0 1 2 3 4 5 6 0 0.5 1 1.5 t (s) y t 1 t N ′ Figure 1.5:
Step responses for some frozen systems obtained every 250 s. The system responses for t 1 = 0 and t N ′ = 9000 s, matching in this case, are highlighted.

Indirect vs direct identification

For the indirect identification we need to determine the discrete-time equivalent of the true system (1.61). Considering the appropriate zero-order hold assumption for the input, in this case, it is given by

x(t k ) + α o 1 (t k )x(t k-1 ) + α o 2 (t k )x(t k-2 ) = β o 1 (t k )u(t k-1 ) + β o 2 (t k )u(t k-2 ) (1.65a) y(t k ) = x(t k ) + e(t k ) (1.65b)
For the model set we assume a DT LTV ARX representation, i.e.

y(t k ) + α 1 (t k )x(t k-1 ) + α 2 (t k )y(t k-2 ) = β 1 (t k )u(t k-1 ) + β 2 (t k )u(t k-2 ) + e(t k ) (1.66)
The model is estimated with the MATLAB routine recursiveARX with a KF as adaptation mechanism. That corresponds to the KF algorithm (1.43). By trial and error, we define the normalized covariance matrix Q n (see (1.46)), which is given by

Q n = diag([10 -3 10 -3 10 -5 10 -5 )
The true system is an OE representation, which is a common description of many physical processes. However, for the model set we have chosen an ARX representation which corresponds to a linear regression. The advantage of the linear regression is that it can be estimated using RLS, which is a simple and computationally efficient algorithm. However, there is a mismatch between the noise of the true system and the model, which would lead to worse estimates. To compensate this effect, data prefiltering will be considered. In the direct approach we will use RLSSVF, with an SVF (see (1.19)) of order n a = 2 and cut-off frequency λ svf = 3.2 (rad/s); note that for λ svf we have conveniently chosen a value slightly larger than the maximum system bandwidth (see Figure 1.7). In the indirect method, for data prefiltering, we consider the same SVF discretized with the bilinear transform. The DT estimates are shown in Figure 1.8; notice that the estimates for α i (t k ) are better than the estimates for β i (t k ). At every step, the conversion from DT to CT parameters is carried out using the appropriate zero-order hold assumption. The transformed CT estimates from the indirect method are presented in Figure 1.9(a).

In the direct approach we assume a CT LTV ARX representation for the model set, i.e. y

(2) Again, there is a mismatch between the noise of the true system and the model set. However in this case, there is no need of data prefiltering since the direct method considers a prefiltering operation in order to obtain the prefiltered time-derivatives. The model is estimated with RLSSVF-KF (1.59). We need to select the normalized covariance matrix Q n . Since the parameters vary at a similar rate, according to Remark (1.3), we can simplified the problem and use Q n = σ 2 w • I 3 ; by trial and error we choose Q n = 10 -4 • I 3 . In Figure 1.9(b), the RLSSVF estimates are plotted.

f (t k ) + a 1 (t k )y (1) f (t k ) + a 2 (t k )y f (t k ) = b 0 (t k )u(t k ) + e(t k ) (1.67)
Both the indirect and direct methods are able to track the CT parameters, even when the mismatch between the noise of the true system and the model introduces errors in the estimates. However, in the indirect approach we estimate first 4 DT parameters while in the direct approach 3. Therefore, in the direct approach Q n is 1 dimension smaller, which means that it is easier to choose than in the indirect approach.

In Section 2.3.5, this example is revisited considering the right OE model and IV based methods. As we can expect, that will lead to better results.

Recursive identification in practice

Applying on-line identification in real-life problems is challenging, because practical and implementation issues have to be considered. In [START_REF] Fridholm | Robust recursive impedance estimation for automotive lithium-ion batteries[END_REF] for instance, some issues like noise effects, model errors and numerical aspects have been addressed in the context of automotive lithium-ion batteries. Literature discussing and proposing solutions to these topics are e.g. [START_REF] Niedźwiecki | Identification of Time-varying Processes[END_REF] and [START_REF] Åstrom | Adaptive Control[END_REF].

A very important practical issue in recursive estimation is the poor excitation situation, which has received a lot of attention in the last decades [START_REF] Cao | A novel recursive algorithm for directional forgetting[END_REF], 2004[START_REF] Stenlund | Avoiding windup in recursive parameter estimation[END_REF][START_REF] Evestedt | Windup properties of recursive parameter estimation algorithms in acoustic echo cancellation[END_REF]. Input richness is required for the stability and convergence of the algorithms. Under poor excitation conditions, the parameter estimates will not be accurate and in the worst case scenario they will blow up. This phenomenon is called estimator windup, and strategies to solve the problem will be presented.

Another practical issue is the choice of hyperparameters of the identification algorithms. For the Kalman filter, the hyperparameters are the covariance matrices that can be estimated using different methods like Maximum Likelihood and Expectation Maximization [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF][START_REF] Bavdekar | Identification of process and measurement noise covariance for state and parameter estimation using extended Kalman filter[END_REF]. In this thesis, the former is adapted to the direct CT identification scenario.

Applications in recursive identification

Regarding implementation issues, we will address the numerical aspects of the estimation algorithms. In batch estimation, standard tools to improve the conditioning of the least squares problem are the singular value decomposition and QR-factorization [START_REF] Ljung | System Identification. Theory for the User[END_REF]. In recursive estimation, which is typically implemented in microcontrollers, this issue may be more critical due to roundoff errors introduced as a consequence of poor machine precision. Several solutions, like square root filtering (see e.g. [START_REF] Grewal | Kalman Filtering[END_REF]), have been proposed to deal with this problem. The use of such methods will be discussed in the particular context of CT identification. Additionally other implementation issues will be tackled.

Applications in recursive identification

Applications in which recursive schemes are used can be divided in the following categories according to [Gustafsson, 2000, p. 4]:

• Surveillance and parameter tracking. For example, monitoring and adaptive control are part of this category.

• State estimation. Observer-based controller, like Linear Quadratic Gaussian (LQG) controllers belongs to this category.

• Fault diagnosis. The use of detection together with isolation techniques, usually called diagnosis, is part of this area.

A precise definition of some of the terms mentioned above, like surveillance and monitoring, that are related to fault detection, can be found in [Gustafsson, 2000, p. 475].

It is important to keep in mind, that the evaluation criteria of the methods depends on the application [START_REF] Gustafsson | Adaptive Filtering and Change Detection[END_REF]. Direct CT identification methods proposed in this thesis are suitable for monitoring physical parameters. The goal in such an application is to have small tracking errors, while in other it might be more important to have small prediction errors.

There are many successful applications in which recursive estimation methods have been used. In control systems for instance, some applications using an ABB adaptive controller are presented in [START_REF] Åstrom | Adaptive Control[END_REF]. In adaptive regulation, where the aim is to attenuate time-varying disturbances, some applications are given in [START_REF] Landau | Adaptive regulation -Rejection of unknown multiple narrow band disturbances (a review on algorithms and applications)[END_REF]. Regarding fault detection/diagnosis, examples of applications are presented in [START_REF] Isermann | Fault diagnosis of machines via parameter estimation and knowledge processing -Tutorial paper[END_REF][START_REF] Sbarbaro | Adaptive soft-sensors for on-line particle size estimation in wet grinding circuits[END_REF].

The application that has motivated our research is an electronic throttle control (ETC), which is a system used in automobiles to regulate air flow in the internal combustion engine. These devices have been intensively studied in the last two decades (see [START_REF] Pavković | Adaptive control of automotive electronic throttle[END_REF][START_REF] Bernardo | Synthesis and experimental validation of the novel LQ-NEMCSI adaptive strategy on an electronic throttle valve[END_REF][START_REF] Zhang | LPV Modeling and mixed constrained H 2 /H ∞ control of an electronic throttle[END_REF] and references therein). In simple terms, the ETC works as follows: the airflow is regulated by a throttle valve which is actuated by a DC servomotor that is regulated by a controller. This system is time-varying because of aging effects and temperature variations. In Chapter 4, our aim is to use the developed techniques to recursively estimate the physical parameters of a throttle valve model, that could be useful for monitoring, an adaptive controller or a fault diagnosis mechanism. Since we are interested in the physical parameters of the system, we focus on direct CT model identification.

Organization of this thesis and contributions

In CT LTI system identification, IV based approaches are an efficient tool to reduce bias in the estimates. In particular, the SRIVC method can deliver optimal (unbiased and minimum variance) estimates. In CT LTV system identification, IV allows us to obtain more accurate models. In Chapter 2, the goal is to extend IV based methods to the estimation of LTV systems operating both in open loop and closed loop. Through numerical examples we show the overall good performance of these algorithms.

In practice, recursive identification of DT time-varying systems is challenging, because several practical aspects and implementation issues have to be considered. In the recursive algorithms for CT LTV systems that we proposed in Chapter 2, similar practical aspects and implementation issues appear. The goal of Chapter 3 is to present suitable solutions to some of them. Regarding practical aspects, we address the problem of estimating hyperparameters of the presented algorithms. Another important aspect is the problem of poor excitation. Implementation issues which lead to improvement in numerical aspects are also discussed.

Chapter 4 is dedicated to two applications. The first one is a real-life electronic circuit which has been recently proposed as a benchmark for both LPV and LTV model identification [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF]. The second one is the throttle valve where it is assumed that the model parameters vary due to the temperature. Since there are nonlinearities in the system, the identification method proposed in Chapter 2 has to be adapted to this particular case.

Publications

Instrumental variable based methods

Introduction

Instrumental variable (IV) methods, studied in statistics, econometrics and engineering, have a long history. One of the first publications in engineering is [START_REF] Wong | Identification of linear discrete time systems using the instrumental variable method[END_REF]. The literature on this technique is extensive (see e.g the books [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF][START_REF] Young | Recursive Estimation and Time-Series Analysis: An Introduction[END_REF][START_REF] Ljung | System Identification. Theory for the User[END_REF] and the references therein). Recent studies are for instance [START_REF] Gilson | Instrumental variable methods for closed-loop system identification[END_REF][START_REF] Laurain | Direct identification of continuous-time linear parameter-varying input/output models[END_REF][START_REF] Tóth | Instrumental variable scheme for closed-loop LPV model identification[END_REF]. The popularity of IV methods lies in the simplicity of some of its variants together with superior statistical properties in comparison with the least squares.

In CT LTI system identification, it is well-known that LSSVF, although simple, is unsatisfactory because the resultant parameter estimates are asymptotically biased. Instrumental variable methods are suitable approaches to deal with bias reduction. Two approaches, mentioned briefly in Chapter 1, are:

• Instrumental variable state-variable filter (IVSVF).

• Simplified refined instrumental variable method of estimation for continuous-time models (SRIVC).

In CT LTV systems, IV based methods allow us to obtain more accurate models. In this chapter, the methods IVSVF and SRIVC are extended to LTV systems, both for open loop and closed loop identification. The contents of this chapter consists basically in two parts; the first one dedicated to open-loop identification and the second to closed-loop identification. The former is partly based on [START_REF] Padilla | Real-time identification of continuoustime linear time-varying systems[END_REF], and it is organized as follows: In Section 2.2, we address the identification of linear time-invariant (LTI) systems. First we recall the optimal IV solution and then the IVSVF and SRIVC methods. Afterwards we present the recursive algorithms for the LTI case. Section 2.3 is dedicated to the identification of linear-time varying models. After formulating the problem, we discuss how to handle the problem of the time-derivatives. Then, the recursive IV based methods for the LTV case are developed. In Section 2.3.4 we discuss measures to assess the estimator performance. We finish the first part with two numerical examples.

The second part, dedicated to the closed-loop case, is based on [START_REF] Padilla | Recursive online IV method for identification of continuous-time slowly time-varying models in closed loop[END_REF]] and it is organized as follows: In Section 2.4.1, the identification problem is formulated. Then, in Section 2.4.2, the optimal closed loop IV solution is revisited. In Section 2.4.3, the closed-loop IV identification methods for LTI models are reviewed. Afterwards, the recursive algorithms for the LTI case are developed. The proposed closed-loop IV based identification methods for linear CT slowly time-varying systems are given in Section 2.4.5 and afterwards a numerical example is presented in Section 2.4.6.

Open-loop identification of LTI models

Optimal IV estimation

Our goal is to find a model for the CT LTI OE system (1.1). In Chapter 1, we propose to represent S with a CT LTI ARX model and estimate it using LSSVF. However, the LSSVF estimate is biased because ϕ f (t k ) and v f (t k ) are correlated [START_REF] Garnier | Direct continuous-time approaches to system identification. Overview and benefits for practical applications[END_REF], since both terms are correlated with e(t k ).

In this thesis, we propose to solve this problem with IV techniques; arguments to favor this choice have been given in Section 1.2.2. Let us recall the optimal IV method where the estimate is defined by θ = arg min

θ∈R n θ 1 N N k=1 F (p)ζ(t k )F (p)ϕ T (t k ) θ - N k=1 F (p)ζ(t k )F (p)y (na) (t k ) 2 W (2.1)
where F (p) is a stable prefilter and x 2 W = x T W x, with W a positive definite weighting matrix. ζ(t k ) is the instrument vector that can be computed in different ways. If S ∈ M, the estimates (2.1) are consistent under the following conditions [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]

1 : C1. Ē{F (p)ζ(t k )F (p)ϕ T (t k )} is full column rank. C2. Ē{F (p)ζ(t k )F (p)v o (t k )} = 0. with v o (t k ) = A o (p)e o (t k ).
Optimal estimates, i.e. unbiased and minimum variance estimates, are obtained if the following additional conditions are satisfied (see [START_REF] Tóth | Instrumental variable scheme for closed-loop LPV model identification[END_REF] and references therein)

C3. W = I C4. F (p) = 1 Ao(p)
C5. The instrument vector is computed using the auxiliary model

ζ(t k ) = -x (na-1) (t k ) . . . -x(t k ) u (n b ) (t k ) . . . u(t k ) T (2.2) where x(t) = B o (p) A o (p) u(t) (2.3)
It is important to notice that the optimal solution requires knowledge of the true system, i.e. the optimal instrument (2.2) and optimal auxiliary model (2.3) are built assuming that the true system is known.

1 The notation Ē[.] = limN→∞ 1 N N k=1 E[.
] is adopted from the prediction error framework of [START_REF] Ljung | System Identification. Theory for the User[END_REF].

Off-line IV algorithms

A basic suboptimal solution to (2.1) is the IVSVF method. This is a two-step approach with F (p) defined in (1.19). In the first step, LSSVF estimates θ1 are computed from (1.32). In the second step, the estimates are obtained using the instrumental variable estimator, i.e.

θ2 = N k=1 ζ f (t k , θ1 )ϕ f T (t k ) -1 • N k=1 ζ f (t k , θ1 )y (na) f (t k ) (2.4)
where y

(na) f (t k ) and ϕ f (t k ) are given in (1.25). The filtered instrument ζ f (t k , θ1 ) is defined through the instrument ζ(t k , θ1 ) = -x (na-1) (t k ) . . . -x(t k ) u (n b ) (t k ) . . . u(t k ) T (2.5) with x(t k )
given by the auxiliary model

x(t k ) = B(p, θ1 ) A(p, θ1 ) u(t k ) (2.6) Then, ζ f (t k , θ1 ) = F (p)ζ(t k , θ1 ) = -x (na-1) f (t k ) . . . -x f (t k ) u (n b ) f (t k ) . . . u f (t k ) T (2.7)
An optimal IV method for CT OE models is the SRIVC method. This is an iterative method, where both the prefilter and the instruments are updated in each iteration based on the parameter estimates obtained at the previous iteration. The SRIVC estimate at the ith iteration is given by θi

= N k=1 ζ f (t k , θi-1 )ϕ f T (t k , θi-1 ) -1 • N k=1 ζ f (t k , θi-1 )y (na) f (t k , θi-1 ) (2.8)
In this case, the filtered linear regression (1.24) is obtained using the following adaptive prefilter

F (p, θi-1 ) = 1 Â(p, θi-1 )
(2.9) Therefore,

y (na) f (t k , θi-1 ) = p na F (p, θi-1 )y(t k ) (2.10a) ϕ T f (t k , θi-1 ) = F (p, θi-1 )ϕ T (t k ) (2.10b) The filtered instrument ζ f (t k , θi-1 ) is defined through the instrument ζ(t k , θi-1 ) = -x (na-1) (t k ) . . . -x(t k ) u (n b ) (t k ) . . . u(t k ) T (2.11)
with x(t k ) defined by the auxiliary model

x(t k ) = B(p, θi-1 ) Â(p, θi-1 ) u(t k ) (2.12)
Figure 2.1: Summary of linear filter methods.

Then,

ζ f (t k , θi-1 ) = F (p, θi-1 )ζ(t k , θi-1 ) = -x (na-1) f (t k ) . . . -x f (t k ) u (n b ) f (t k ) . . . u f (t k ) T (2.13)
The three linear filter methods presented so far are summarized in Figure 2.1. We can see for each approach the corresponding prefilter and auxiliary model.

Remark 2.1

In a more general case, the measurement noise in S might be colored. A reliable solution is then to use the RIVC method [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF]. Nevertheless, the SRIVC algorithm has the advantage of its much simpler implementation in this recursive estimation context and it performs optimally if the additive noise is white, i.e. estimates are consistent and asymptotically efficient. Moreover, if the noise is colored, with rational spectral density, then the estimates remain consistent, thanks to the inherent instrumental variable mechanism; and although they do not have minimum variance, experience has shown that the estimates are often relatively efficient ( i.e. they have low but not minimum variance [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF][START_REF] Garnier | Direct continuous-time approaches to system identification. Overview and benefits for practical applications[END_REF]). Taking into account the trade-off between accuracy and implementation issues, the recursive version of both IVSVF and SRIVC are considered in this work.

Recursive instrumental variable methods

The recursive version of IVSVF (RIVSVF) is the following algorithm, which is based on the methodology used in the DT case, i.e. in this case considering the recursive instrumental variable method (see e.g. [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF]):

θ(t k ) = θ(t k-1 ) + L(t k )ε(t k ) (2.14a) ε(t k ) = y (na) f (t k ) -ϕ f T (t k ) θ(t k-1 ) (2.14b) L(t k ) = P (t k-1 )ζ f (t k ) 1 + ϕ f T (t k )P (t k-1 )ζ f (t k ) (2.14c) P (t k ) = P (t k-1 ) -L(t k )ϕ f T (t k )P (t k-1 ) (2.14d)
where y

(na) f (t k ) and ϕ f (t k ) are defined in (1.25). However, note that the CT filtering operations involved in (1.25) in the off-line case are different than in the recursive case. This implementation issue is discussed later in Section 3.1.3. The filtered instrument

ζ f (t k ) is defined through the instrument ζ(t k , θ(t k-1 )) = -x (na-1) (t k ) . . . -x(t k ) u (n b ) (t k ) . . . u(t k ) T (2.15)
with x(t k ) given by the auxiliary model

Â(p, θ(t k-1 ))x(t k ) = B(p, θ(t k-1 ))u(t k ) (2.16) Then, ζ f (t k ) = F (p)ζ(t k , θ(t k-1 )) = -x (na-1) f (t k ) . . . -x f (t k ) u (n b ) f (t k ) . . . u f (t k ) T (2.17)
For the sake of simplicity, hereafter and whenever necessary, the parameter dependency will be sometimes omitted, as it is the case for ζ f (t k ) in (2.17), which should be

ζ f (t k , θ(t k-1 )
). Contrary to the off-line case, where the dependency is w.r.t. an estimate of a previous step, here it is w.r.t. an estimate of a previous recursion.

For the recursive version of SRIVC (RSRIVC), the algorithm is given by (2.14), but with

y (na) f (t k ) = p na F (p, θ(t k-1 ))y(t k ) (2.18a) ϕ f T (t k ) = F (p, θ(t k-1 ))ϕ T (t k ) (2.18b)
where the adaptive prefilter is

F (p, θ(t k-1 )) = 1 Â(p, θ(t k-1 )) (2.19) From ζ(t k ) given in (2.15), we define the filtered instrument ζ f (t k ) by ζ f (t k ) = F (p, θ(t k-1 ))ζ(t k , θ(t k-1 )) = -x (na-1) f (t k ) . . . -x f (t k ) u (n b ) f (t k ) . . . u f (t k ) T (2.20)
It should be noted that in this case, not only x(t k ), but also the prefilter (2.19) depend on the estimate θ(t k-1 ), meaning actually that all prefiltered variables depend on θ(t k-1 ). Notice also that in off-line identification with SRIVC, the auxiliary model and prefilter depend on an estimate of a previous iteration, whereas in the recursive case they depend on an estimate of a previous recursion.

Both RIVSVF and RSRIVC are initialized with RLSSVF. The transition from the latter to the IV based methods can be manually set by the user at a time instant that we denote t s,iv . Techniques that allow us to perform this transition automatically will be discussed in Chapter 3.

Open-loop identification of LTV models

In this section, the aim is to extend these IV approaches presented previously to realtime estimation of linear time-varying systems.

Problem formulation

Let us consider a CT LTV OE system with input u and output y. The data generating system is given by

S A o (p, t)x(t) = B o (p, t)u(t) y(t k ) = x(t k ) + e o (t k ) (1.49)
where e o (t k ) is a zero-mean DT white noise sequence. We assume that the system (1.49) can be represented by the model set

M A(p, t, θ)x(t) = B(p, t, θ)u(t) y(t k ) = x(t k ) + e(t k ) (2.21)
A(p, t, θ) and B(p, t, θ) are the following polynomials with time variable parameters:

B(p, t, θ) = b 0 (t)p n b + b 1 (t)p n b -1 + . . . + b n b (t) A(p, t, θ) = p na + a 1 (t)p na-1 + . . . + a na (t)
where n a ≥ n b and e(t k ) is a zero-mean DT white noise. The time-varying parameters can be gathered in the parameter vector θ(t), i.e.

θ(t) = a 1 (t) . . . a na (t) b 0 (t) . . . b n b (t) T (2.22)
The following assumptions are supposed to be satisfied: A1. The system and the model set have CT LTV OE representation, with polynomial degrees n a and n b a priori known.

A2. The true parameter vector θ o (t) is slowly time-varying in the sense that in a local window with time interval s ∈ [t i , t f ], it can be locally approximated by a constant parameter vector θ c , i.e.

θo (t) ≤ ǫ θ such that θ o (s) -θ c ≤ ǫ ∆θ (2.23)
where ǫ θ and ǫ ∆θ are small numbers.

A3. The input u is persistently exciting.

Then, the identification problem is to recursively estimate the time-varying parameters that characterize the model structure given by (2.21), based on sequential samples of input and output data

Z N ′ = {u(t k ); y(t k )} N ′ k=1
, where N ′ is the number of samples which increases by one with every recursion.

Handling of the time-derivative issue

As it was briefly mentioned in Section 1.3.3, multiplication with the differentiation operator p is not commutative with the time-dependent variables, i.e.

pB(p, t, θ)u(t) = B(p, t, θ)pu(t)

Thus, it would appear at first that we are no longer able to filter the LTV model (2.21) without introducing errors and obtain a linear regression, as in the LTI case (see (1.24)). One solution would be to model the time-varying parameters in a deterministic way such that the new parameters are constant or nearly constant, as in [START_REF] Jiang | A recursive identification method for continuous time-varying linear systems[END_REF]. For on-line estimation, the difficulty of such an approach lies in finding suitable deterministic models. Moreover, extra parameters are introduced and, as a consequence, larger variance of the estimates are obtained.

Alternatively, we consider Assumption (A2) in Section 2.3.1, which states that the true time-varying parameters are slowly varying in the sense that the system can be locally approximated by a constant parameter vector, i.e. an LTI model. Recursive estimation algorithms capable of tracking time-varying parameters have a finite memory, meaning that the data used for identification is a sliding window in which the approximation of the LTI system is valid. Since an LTI representation does not suffer from the non-commutativity problem, then the non-commutativity aspect is negligible in the estimation of slowly LTV systems. Thus, as in the LTI case, prefilter time-derivatives can be obtained by applying a low-pass filter to the model (2.21). Analogous to the case for CT LTV ARX models (see (1.53)), applying a filter F (p) to (2.21) we obtain

y (na) f (t k ) = ϕ f T (t k )θ(t k ) + v f (t k ) (2.24)
where y

(na) f (t k ) and ϕ f (t k ) are defined in (1.25). The difference between (1.53) and (2.24) is that in the latter the filtered noise is given by

v f (t k ) = F (p)A(p, t, θ)e(t k ) (2.25) Remark 2.2
In [START_REF] Dimogianopoulos | Adaptive control for linear slowly time-varying systems using direct least-squares estimation[END_REF], the same reasoning is done to compute prefilter time-derivatives and estimate a CT model. Moreover, the authors provided Lemma 1 which guarantees the existence of a local LTI representation under the conditions of a bounded change ratio of the true parameters and bounded measurement noise. Lemma 1 has been derived from a exponentially weighted least squares in continuous-time in order to have an explicit dependency on the change ratio of the true parameters. However, Lemma 1 can be also expressed in DT using approximations of the derivatives and integrals involved in the formulation.

In practice, a simple way to validate Assumption (A2) would be to compare the cost functions

V [ θ(s)] t i ≤s≤t f = 1 N f k=i [y (na) f (t k ) -ϕ f T (t k ) θ(t k )] 2 (2.26) with V [θ c ],
where θ c ∈ θ(s). Then, Assumption (A2) is valid, if

V [θ c ] V [ θ(s)] t i ≤s≤t f ≈ 1 (2.27)
Remark 2.3 In tracking analysis it can be distinguished between transient mode and tracking mode [START_REF] Rutström | Simplified Wiener LMS Tracking with Automatic Tuning of the Step-Size[END_REF]. The former mode starts when the parameters jump due to either an inaccurate initialization of the recursive algorithm or an abrupt change in the system. After the transient mode, we arrive at the tracking mode which is a steadystate condition. In recursive estimation of LTI models, the estimated parameters vary abruptly in the transient mode, and then they vary slowly in the tracking mode. Thus, in the tracking mode, the non-commutativity can be neglected and an SVF can be applied to circumvent the time-derivative issue. However, in the transient mode this argument is not valid and errors in the computed estimates are introduced.

Recursive instrumental variable methods

For the reasons exposed in Remark 2.1, let us now adapt IVSVF and SRIVC to LTV systems.

To cope with the parameter variation we assume that they can be represented by a random walk model. Then, the full model becomes: The parameters of model (2.28) can be estimated based on the methodology used in the DT case, i.e. in this case considering the instrumental variable based Kalman filter (see [Young, 2011, p. 298]). The algorithm is given by Prediction step:

M        θ(t k ) = θ(t k-1 ) + w(t k ) A(p, t, θ)x(t) = B(p, t, θ)u(t) y(t k ) = x(t k ) + e(t k ) (2.28)
θ(t k |t k-1 ) = θ(t k-1 )
(2.29a)

P (t k |t k-1 ) = P (t k-1 ) + Q n (2.29b)
Correction step:

θ(t k ) = θ(t k |t k-1 ) + L(t k )ε(t k ) (2.29c) ε(t k ) = y (na) f (t k ) -ϕ f T (t k ) θ(t k |t k-1 ) (2.29d) L(t k ) = P (t k |t k-1 )ζ f (t k ) 1 + ϕ f T (t k )P (t k |t k-1 )ζ f (t k )
(2.29e)

P (t k ) = P (t k |t k-1 ) -L(t k )ϕ f T (t k )P (t k |t k-1 ) (2.29f)
where y Another way to deal with the estimation of time-varying parameters is the forgetting factor approach. Based on the methodology used in the DT case, the forgetting factor based RIVSVF/RSRIVC algorithm is defined by

(na) f (t k ), ϕ f (t k
θ(t k ) = θ(t k-1 ) + L(t k )ε(t k ) (2.30a) ε(t k ) = y (na) f (t k ) -ϕ f T (t k ) θ(t k-1 ) (2.30b) L(t k ) = P (t k-1 )ζ f (t k ) λ + ϕ f T (t k )P (t k-1 )ζ f (t k )
(2.30c)

P (t k ) = 1 λ P (t k-1 ) -L(t k )ϕ f T (t k )P (t k-1 ) (2.30d)
It should be however reminded that a drawback of those approaches is that they are not suitable for cases where parameters vary at different rates. Moreover, as we will see in Section 3.6.1, KF based approaches are in general more robust to poor excitation. This motivates the focus on KF based methods in this work.

Both the forgetting factor and Kalman filter based routines are available in the CONtinuous-Time System IDentification (CONTSID) toolbox, which has been recently updated [START_REF] Padilla | Version 7.0 of the CONTSID toolbox[END_REF] and can be freely downloaded2 .

Measures of estimator performance

In practice, the tracking ability of the approaches can be measured only indirectly in terms of the prediction error ε. Nonetheless, in simulation the true parameters are available, which allows us to assess the performance by means of the tracking error

θ(t k ) = θ o (t k ) -θ(t k ) (2.31)
For N samples, the mean squared value of the parameter tracking error is given by

MSE[ θ(t k )] = 1 N N k=1 θ o (t k ) -θ(t k ) 2 (2.32)
For the ith element of θ o (t k ), denoted by θ o i (t k ), we can also use the relative error defined by θi,r

(t k ) = θ o i (t k ) -θi (t k ) θ o i (t k ) × 100 (2.33)
By adding and subtracting E[ θ(t k )] in (2.31), the tracking error can be split into two parts [Haykin, 2014, p. 524], the lag error θl (t k ) and the noise error θn (t k ), i.e.

θ(t k ) = θ o (t k ) -E[ θ(t k )] θl (t k ) + E[ θ(t k )] -θ(t k ) θn(tk)
(2.34) These errors cannot be minimized simultaneously due to the well-known trade-off of recursive algorithms between tracking capability and noise sensitivity [START_REF] Ljung | Adaptation and tracking in system identification -a survey[END_REF]. For tracking problems, the lag error can be compared to the standard concept of bias in LTI systems [Rutström, 2005, p. 15]. The presence of the lag error is testimony to the non-stationary nature of the environment [START_REF] Haykin | Adaptive Filter Theory[END_REF], meaning that in the case of time-varying parameters, the lag error is larger than in the time-invariant case. Indeed, in the latter, the lag error could be zero.

The lag error and noise error could be computed by a Monte Carlo simulation using a given input and different measurement noise realizations. Then, the expectation E[ θ(t k )] will be replaced by the sample mean θ(t k ) taken over the total number of simulations.

Such an analysis is done in Section 2.3.5. It is then convenient to define two additional measures, the mean squared value of the lag error defined by

MSE l [ θ(t k )] = 1 N N k=1 θ o (t k ) -θ(t k )] 2 (2.35)
and the mean squared value of the noise error given by

MSE i n [ θ(t k )] = 1 N N k=1 θ(t k ) -θ(t k ) 2 (2.36)
where the superscript i means that the computation is done w.r.t. the estimate θ(t k ) from the ith simulation which can be arbitrarily chosen.

Remark 2.4 In the identification of LTI models, IV methods are applied to reduce bias. Analogously, in the LTV case, the aim is to reduce the lag error which come at the expenses of increasing slightly the noise error. That is the well-known trade-off of recursive algorithms between tracking ability and noise sensitivity.

Numerical example -Similar parameter variations

Example 1.3.4 is revisited here to illustrate the performance of the proposed IV based Kalman filter methods. Additionally, the lag error and noise error defined in Section 2.3.4 are analyzed. For all the simulations, we use the same value for the hyperparameters Q n and λ svf . The former is a diagonal matrix, and since the parameters are varying at a similar rate we consider

Q n = 10 -4 I 3
where I 3 is the identity matrix of dimension 3. For the SVF, λ svf = 3.2 rad/s is chosen, i.e. a convenient value slightly larger than the maximum bandwidth (see Figure 1.7). The IV methods, which are initialized with RLSSVF, start to operate at time instant t s,iv = 100 s. The value for t s,iv is chosen based on the convergence of both the output error and the parameter tracking error.

We present next some results for a single experiment run out of the 100 that are considered later in a Monte Carlo simulation analysis.

Single experiment analysis

In the LTI case, it is known that the RLSSVF estimates are always biased due to the measurement noise. Even if the bias cannot be fully removed, it can be reduced by a proper choice of the cut off frequency λ svf . In the LTV case, this is more difficult since the system bandwidth is varying, while the SVF bandwidth is constant. Instrumental variable methods, like RSRIVC, can be used to cope with the noise and reduce the lag error, as explained in Remark 2.4. The RLSSVF and RSRIVC estimates are compared with the true parameter in Figure 2.3, from where we can see the improvement obtained by using the IV approach. The mean squared value of the parameter tracking error is compared in Table 2.1, where we can see that the performance of RIVSVF is similar to RSRIVC, the latter being slightly better as we could expect, due to the adaptive prefiltering. Note that the value for λ svf has been conveniently chosen; an unsuitable cut-off frequency would produce worse results.

Monte Carlo simulation analysis

The performance of these algorithms is further assessed through a Monte Carlo simulation analysis with 100 experiments. As a first performance measure, we take the relative error θi,r (2.33). Averaging the relative errors over the 100 experiments, we obtain the results presented in Figure 2.4 for RLSSVF and RSRIVC. The improvement is clear from the plot.

For LTV model identification, the purpose of using IV methods is to reduce the lag error (2.34) as we can see from Figure 2.5, where its absolute value is shown for RLSSVF and RSRIVC. The absolute value of the noise error is plotted in Figure 2.6 only for the last 1000 s, to be able to distinguish the RLSSVF result from the RSRIVC result. The reduction of the lag error implies an increase of the noise error, although in terms of error magnitudes, the reduction of the former is greater than the increase in the latter. This is confirmed from Table 2.2, where the mean squared of the lag error (2.35) and mean squared of the noise error (2.36) w.r.t. experiment 1, are presented. Here we can see additionally the results for RIVSVF, whose MSE l [ θ(t k )] is slightly worse than the one for RSRIVC as expected.

Numerical example -Different parameter variations

To illustrate further the performance of the proposed algorithms, let us consider a particular LTV system whose bandwidth variation in time is significant. The data generating system is now given by the following second order CT LTV OE system. The sampling time is set to 0.01 s and the total simulation time is 1000 s. The input is a PRBS and the DT measurement noise is a zero-mean, Gaussian noise with constant variance 0.1. Part of the input-output data is shown in Figure 2.7.

S p 2 + a o 1 (t)p + a o 2 (t) x(t) = b o 0 (t)u(t) y(t k ) = x(t k ) + e(t k ) (2.37)
In this example, as a consequence of the time-varying parameters, the DC gain is decreasing towards half of the simulation time; and since the noise variance is kept constant, the signal-to-noise ratio (SNR) is decreasing around half of the simulation time. For some frozen systems obtained every 100 s, the Bode diagrams are plotted in Figure 2.8, where we can see how the DC gain is varying. The step responses and bandwidths of these frozen systems are shown in Figures 2.9 and 2.10, respectively. Notice that the ratio between the maximum and minimum bandwidths is nearly 10, i.e. the bandwidth variation is relatively large over the total simulation time, as we mentioned.

For all the estimations in this example, the same normalized covariance matrix Q n is used. Q n has been obtained by a numerical search; it is set as a diagonal matrix with diagonal elements equal to 10 -5 , 10 -4 and 10 -10 corresponding to a o 1 , a o 2 and b o 0 , respectively. The value corresponding to b o 0 is 10 -10 , because it is assumed known that this parameter is constant. For the SVF, λ svf = 16 rad/s is chosen, i.e. a value slightly greater than the maximum bandwidth. The IV based methods are initialized with RLSSVF. Based on the convergence of RLSSVF, the start of the IV methods t s,iv is set to 10 s. In addition, for the IV methods we consider the filter This filter is used to feed back the algorithms with the estimated parameters (see more details in Section 3.2.3). We present next some results for a single experiment run out of the 100 that are considered later in a Monte Carlo simulation analysis.

M θ (q -1 ) = 0.04762 + 0.04762q -1 1 -0.9048q -1

Single experiment analysis

Firstly, the RLSSVF results are analyzed. The RLSSVF estimates are displayed in Figure 2.11 for two situations: the first a noise-free (or deterministic) case, and the second using noisy data. In the former case, since there is no noise, RLSSVF performs very well and it is difficult to distinguish between the true and estimated parameters. However, in the noisy case, RLSSVF does not provide reasonable parameter estimates, as it can be observed in Figure 2.11(b). Note that the relative errors of the three estimates are similar but this is difficult to see because the limits for the y-axes are different. Similar results have been obtained for other values of λ svf . Notice that, in this example, the variance of the tracking error increases slightly over the central portion of the plot that can be seen in Figure 2.11. This is because the SNR is decreasing in time while the system bandwidth reduces in time (which would also require λ svf to be lower).

To cope with the measurement noise and reduce the lag error, RIVSVF and RSRIVC are used. The improvement in performance obtained by using these IV based methods is clear from Figure 2.12. We can also see that the RSRIVC algorithm gives better results than RIVSVF at the end of the simulation time as expected, because the RSRIVC prefilter is adapted over time. Naturally, the changes of the variance of the tracking error observed in the RLSSVF case, also occur here.

Monte Carlo simulation analysis

A more thorough analysis is done by performing a Monte Carlo simulation with 100 experiments. For RLSSVF, the average value of the relative error over the 100 experiments is shown in Figure 2.13. The same is done for the IV methods in Figure 2.14, where we can see clearly the improvements in the results. Notice also that RSRIVC is outperforming RIVSVF for the last 200 seconds. The worse results for RIVSVF are expected since in this part of the simulation, the chosen value for λ svf (16 rad/s) is significantly larger than the system bandwidth (see Figure 2.10).

Closed-loop identification

In real-life, experiments for model identification are often done in closed loop because there are economic or safety constraints, or the system is unstable. The issue with closedloop identification is that the feedback mechanism introduces correlation between the disturbances and the input signal, yielding biased estimates if no special strategy is used. For LTI systems, instrumental variable (IV) techniques have been proposed to solve this problem both for DT models [Gilson andVan den Hof, 2005, Gilson et al., 2011] and CT models [START_REF] Gilson | Instrumental variable methods for closed-loop continuous-time model identification[END_REF]; [see also, Chapter 9 in [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF]. IV have been also extended to the identification of DT linear parameter-varying models in [START_REF] Tóth | Instrumental variable scheme for closed-loop LPV model identification[END_REF]. The recursive identification schemes developed in the previous sections are adapted here for identifying linear CT slowly time-varying systems operating in closed loop.

Problem formulation

Let us consider the closed-loop configuration from Figure 2.15 with a CT LTV OE system S,

S A o (p, t)x(t) = B o (p, t)u(t) y(t k ) = x(t k ) + e o (t k ) (2.38)
as in the open-loop case. G o is the plant and e o (t k ) a zero-mean DT white noise sequence.

In Figure 2.15, the CT controller C can be any nonlinear and/or time-varying controller.

Knowing C, we can compute the input u(t) as follows

u(t k ) = r 1 (t k ) + C(r 2 (t k ) -y(t k )) (2.39)
where C is the operator form of the controller and r 1 (t k ), r 2 (t k ) are external signals. We assume that the system can be represented by the model set M given in (2.21). Suppose that assumptions A1 and A2 defined in Section 2.3.1 for the open-loop case are satisfied. Additional assumptions are: A3. The controller C is known.

A4. The controller C ensures BIBO stability of the closed-loop system defined by (2.38) and (2.39).

A5. The reference signal r(t

k ) = r 1 (t k ) + Cr 2 (t k ) is persistently exciting.
Then, the identification problem is to recursively estimate the time-varying parameters that characterize the model structure given by (2.21), based on the data set

Z N ′ = {r(t k ), u(t k ), y(t k )} N ′ k=1
where N ′ is the number of samples which increases by one with every recursion.

Optimal off-line IV estimation of LTI models

In this section we assume that both the plant and the controller are LTI systems, i.e.

G o : G o (p) = B o (p) A o (p) (2.40) C : C c (p) = Q c (p) P c (p) (2.41)
with the pairs (A o , B o ) and (P c , Q c ) assumed to be coprime. First we address the off-line estimation problem and afterwards its recursive counterpart.

When the parameters in (2.21) are constant, they could be estimated as in the open-loop case using the LSSVF method (see (1.32)). However, the parameters will be biased because ϕ f (t k ) and v f (t k ) are correlated due to the filter and because v f (t k ) is colored. Additionally, in a closed-loop configuration, ϕ f (t k ) and v f (t k ) are correlated due to the feedback mechanism. A solution is to use the closed-loop IV method given by (2.1) (see [START_REF] Gilson | Instrumental variable methods for closed-loop continuous-time model identification[END_REF]). The only difference w.r.t. to the open-loop case is the computation of the instruments. Then, optimal estimates, i.e. unbiased and minimum variance estimates, can be obtained if the conditions C1-C4 defined in Section 2.2.1 together with the following condition C5 are satisfied [START_REF] Gilson | Instrumental variable methods for closed-loop continuous-time model identification[END_REF] C5. The instrument vector is computed using the auxiliary model from Figure 2.16 as follows

ζ(t k ) = -x (na-1) (t k ) . . . -x(t k ) ů(n b ) (t k ) . . . ů(t k ) T (2.42) where A o (p)x(t k ) = B o (p)ů(t k ) (2.43a) ů(t k ) = r 1 (t k ) + C(r 2 (t k ) -x(t k )) (2.43b)

Off-line IV algorithms of LTI models

Conditions C1-C3 can be readily satisfied [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]]. However C4-C5 require knowledge of the true unknown system, which is the usual dilemma with optimal estimation. To circumvent this problem, different estimation algorithms have been proposed, which vary depending on the choice of the instruments ζ(t k ), the filter F (p) and the model structure of the true system [START_REF] Gilson | Instrumental variable methods for closed-loop continuous-time model identification[END_REF]. In the following we recall the closed-loop version of IVSVF and SRIVC presented in Section 2.2.2. The closed-loop version of IVSVF, denoted by CLIVSVF, consists of two steps. First, LSSVF estimates θ1 are obtained using (1.32). In the second step, given N samples, the CLIVSVF estimates are computed from (see [START_REF] Gilson | Instrumental variable methods for closed-loop continuous-time model identification[END_REF])

θ2 = N k=1 ζ f (t k , θ1 )ϕ f T (t k ) -1 • N k=1 ζ f (t k , θ1 )y (na) f (t k ) (2.4)
The difference w.r.t. the open-loop case stands in the computation of the filtered instrument ζ f (t k , θ1 ), which is defined through the instrument

ζ(t k , θ1 ) = -x(na-1) (t k ) . . . -x(t k ) û(n b ) (t k ) . . . û(t k ) T (2.44)
with x(t k ) and û(t k ) given by the auxiliary model

Â(p, θ1 ) x(t k ) = B(p, θ1 ) û(t k ) (2.45a) û(t k ) = r 1 (t k ) + C(r 2 (t k ) -x(t k )) (2.45b) Then, ζ f (t k , θ1 ) = F (p)ζ(t k , θ1 ) = - x(na-1) f (t k ) . . . -xf (t k ) û(n b ) f (t k ) . . . ûf (t k ) T (2.46)
In [START_REF] Gilson | Instrumental variable methods for closed-loop continuous-time model identification[END_REF], a similar algorithm, called CLIVC2, has been proposed. The difference is that in CLIVC2, a CT ARX model instead of a CT OE model, is considered. The closed-loop version of SRIVC, denoted by CLSRIVC, has been developed in [START_REF] Gilson | Instrumental variable methods for closed-loop continuous-time model identification[END_REF]. Given N samples, the CLSRIVC estimates at the ith iteration is θi

= N k=1 ζ f (t k , θi-1 )ϕ f T (t k , θi-1 ) -1 • N k=1 ζ f (t k , θi-1 )y (na) f (t k , θi-1 ) (2.47)
As in the open-loop case, y (na) f

(t k , θi-1 ) and ϕ f (t k , θi-1 ) are obtained from (2.10) using the adaptive prefilter (2.9). The difference w.r.t. the open-loop case stands in the computation of the filtered instrument ζ f (t k , θi-1 ), which is defined through the instrument

ζ(t k , θi-1 ) = -x(na-1) (t k ) . . . -x(t k ) û(n b ) (t k ) . . . û(t k ) T (2.48)
with x(t k ) and û(t k ) given by the auxiliary model

Â(p, θi-1 ) x(t k ) = B(p, θi-1 ) û(t k ) (2.49a) û(t k ) = r 1 (t k ) + C(r 2 (t k ) -x(t k )) (2.49b) Then, ζ f (t k , θi-1 ) = F (p, θi-1 )ζ(t k , θi-1 ) = - x(na-1) f (t k ) . . . -xf (t k ) û(n b ) f (t k ) . . . ûf (t k ) T (2.50)

Recursive estimation of LTI models

The two closed-loop identification methods presented previously can be adapted for recursive estimation. As in the open-loop scenario, the algorithm is given by (2.14).

The difference is only in how the filtered instrument ζ f (t k ) is computed.

For the recursive version of CLIVSVF, denoted by CLRIVSVF, the filtered instrument is defined through the instrument

ζ(t k , θ(t k-1 )) = -x(na-1) (t k ) . . . -x(t k ) û(n b ) (t k ) . . . û(t k ) T (2.51)
with x(t k ) and û(t k ) given by the auxiliary model

Â(p, θ(t k-1 )) x(t k ) = B(p, θ(t k-1 )) û(t k ) (2.52a) û(t k ) = r 1 (t k ) + C(r 2 (t k ) -x(t k )) (2.52b)
where C is the controller. Then,

ζ f (t k ) = F (p)ζ(t k , θ(t k-1 )) = - x(na-1) f (t k ) . . . -xf (t k ) û(n b ) f (t k ) . . . ûf (t k ) T (2.53)
For the recursive version of CLSRIVC, denoted by CLRSRIVC, the filtered instrument is given by

ζ f (t k ) = F (p, θ(t k-1 ))ζ(t k , θ(t k-1 )) = - x(na-1) f (t k ) . . . -xf (t k ) û(n b ) f (t k ) . . . ûf (t k ) T (2.54) with ζ(t k , θ(t k-1
) defined in (2.51) and the prefilter in (2.19).

Recursive estimation of LTV models

In a closed-loop configuration, tracking time-varying parameters changing at different rates can be done using the Kalman filter based RLSSVF. For the reason exposed above, these estimates would be biased

Numerical example

Three recursive algorithms are evaluated: CLRIVSVF, CLRSRIVC and RLSSVF, which is used to initialized the first two approaches. The example presented in Section 2.3.6 is adapted here to the closed-loop configuration shown in Figure 2.15. Regarding the external signals (see Figure 2.15), r 1 (t) is a PRBS and r 2 (t) = 0. e o (t k ) is a zeromean DT Gaussian noise with constant variance 0.1. Remember that the ratio between the maximum and minimum bandwidths is nearly 10, i.e. the bandwidth variation is relatively large over the total simulation time.

CT filtering operations for the computation of prefiltered time-derivatives and instruments are implemented from their discretized counterparts. The discretized version of the PID controller is given by

C : C c (q -1 ) = k p + k i 1 -q -1 + k d (1 -q -1 )
(2.55) with q -1 the backward shift operator and k p = 1.79, k i = 13.8, k d = 5.83 • 10 -2 . For all the simulations, we use the same value for the hyperparameters Q n and λ svf as before.

We present next some results for a single experiment run out of the 100 that are considered later in a Monte Carlo simulation analysis.

Single experiment analysis

In the LTI open-loop case, it is known that the RLSSVF estimates are always biased due to the noise. Even if the bias cannot be removed, it can be reduced by a proper choice of the cut off frequency λ svf . In the open-loop LTV case, this is more difficult since the system bandwidth is varying, while the SVF bandwidth is constant, as it was illustrated in Section 2.3.6. When the system operates in closed loop there is an additional issue, namely the correlation between the input u(t k ) and the noise e o (t k ) due to the feedback mechanism. In the LTI case, that might result in larger bias on the estimates, while in the LTV case, it results in larger lag errors. In order to illustrate the impact of the noise on RLSSVF in the closed-loop LTV situation, we compare the noise-free case with the noisy case. From Figure 2.17, we can see that RLSSVF is able to track the parameters only in the noise-free situation. It is important to stress that the value used for λ svf is a very good choice since it is slightly higher than the maximum system bandwidth.

To circumvent the parameter tracking problem, the two closed-loop IV approaches introduced in Section 2.4.5 have been used. From Figure 2.18 we can see that both CLRIVSVF and CLRSRIVC are able to track the parameters in this situation where the parameters vary at different rates. The former approach considers a constant prefilter while in the latter the filter is automatically adapted from the previous estimates.

Monte Carlo simulation analysis

In order to see more clearly the difference between these algorithms, a Monte Carlo simulation analysis with 100 experiments is run. For comparison purposes, the relative error (2.33) is used as a measure. Averaging the relative errors of the 100 experiments, we obtain the results presented in Figure 2.19 for RLSSVF and in Figure 2.20 for CLRIVSVF and CLRSRIVC. It can be seen that towards the end of the simulation time, the best results are obtained with CLRSRIVC.

Conclusions

In this chapter, the algorithms RIVSVF and RSRIVC have been developed for open-loop identification. Closed-loop versions have been also presented. The recursive algorithms follow from the equivalent discrete-time versions but are more complex because they require prefiltering for handling the time-derivatives. These approaches are able to track slow variations in the model parameters, based on an assumed stochastic random walk model for the parameter variations. One of the advantage of the random walk is that it is suitable when the parameters vary at different rates. RIVSVF can provide good estimates if prior knowledge exists to define the cut-off frequency of the SVF, meaning that it can be chosen close to the true system bandwidth. However, if the system bandwidth vary too much, it will not be possible to choose a value which is suitable for every moment. With RSRIVC, this is not an issue since it considers an adaptive filter which is updated using previous estimates.

The approaches have been evaluated by means of Monte Carlo simulation analysis that demonstrates the effectiveness of the IV approaches to reduce the tracking error, in comparison with RLSSVF. Chapter 3

Practical aspects and implementation issues

Practical aspects and implementation issues in recursive estimation are of crucial importance, and have been discussed e.g. in the books [START_REF] Åstrom | Adaptive Control[END_REF]] and [START_REF] Niedźwiecki | Identification of Time-varying Processes[END_REF]. In this chapter we investigate these problems for the approaches developed in the previous chapter. The particularity of those methods is basically that they handle CT models and they use IV techniques.

The first part of this chapter is dedicated to the choice of the hyperparameters for the methods developed in Chapter 2. We will focus on RSRIVC, since the analysis for the other approaches is analogous. A diagram of RSRIVC highlighting the practical aspects and implementation issues is presented in Figure 3.1; the issues, denoted by Ii in the figure, are the following: I1. Choice of user hyperparameters, which are:

-λ svf : cut-off frequency of the RLSSVF which is used to initialized RSRIVC.

In Section 3.1, aspects related to the SVF, including the choice of λ svf are discussed.

-t s,iv : time instant when RSRIVC starts to operate. In Section 3.2 aspects related to IV techniques, including the choice of t s,iv , are presented.

-Q n : normalized covariance matrix of the Kalman filter. In Section 3.3, an estimation of Q n using maximum likelihood is proposed.

I2. Digital implementation of the CT filtering operations involved in the prefilter and auxiliary model. This is covered in Section 3.1.3.

I3. Numerical aspects. One way to solve this implementation issue is to use square root filtering. For our methods we will see that a suitable approach is to use signal scaling. This is done in Section 3.5.

I4. Estimator windup problem. One of the assumptions made in the theoretical developments in the previous chapter is that the input u(t k ) is persistently exciting. In practice, there might be intervals where this condition is not met. Solutions for this phenomenon, called estimator windup, will be discussed in Section 3.6.

These and some other issues will be discussed along the chapter. For instance, another practical aspect presented in Section 3.4 is how to handle the unknown initial conditions of the system that we want to identify. We will show that by properly handling these initial conditions, the overshoots that they can produce in the estimates during the initial recursions can be reduced.

Aspects related to the SVF

Normalized vs ordinary SVF

In the direct identification approaches presented in Chapter 2, the SVF allows us to obtain prefiltered time-derivatives. For the SVF, two forms are possible, the ordinary version (OSVF) defined in (1.19) and the normalized version (NSVF)

F nsvf (p) = λ svf p + λ svf na (3.1)
Thus, these forms are linked as follows

F nsvf (p) = λ na svf F osvf (p) (3.2)
For the OSVF, given N ′ samples, the parameter covariance matrix is, up to a scaling factor, defined by

P osvf (t N ′ ) = N ′ k=1 ϕ f (t k )ϕ f T (t k ) -1 (3.3)
with ϕ f (t k ) defined in (1.27). Using the NSVF, the filtered regressor vector is given by

ϕ T f,nsvf (t k ) = F nsvf (p)ϕ T (t k ) = λ na svf ϕ f T (t k ) (3.4)
Then, in this case the parameter covariance matrix is, up to a scaling factor, defined as 3.1. Aspects related to the SVF follows

P nsvf (t N ′ ) = N ′ k=1 ϕ f,nsvf (t k )ϕ T f,nsvf (t k ) -1 = 1 λ 2na svf N ′ k=1 ϕ f (t k )ϕ T f (t k ) -1
(3.5) Thus, from (3.3) and (3.5) we obtain [START_REF] Garnier | Influence of the initial covariance matrix on recursive LS estimation of continuous models via generalised Poisson moment functionals[END_REF]]

P osvf (t N ′ ) = λ 2na svf P nsvf (t N ′ ) (3.6)
To use RLSSVF, either with OSVF or NSVF, the initial estimates ( θ(t 0 ), P (t 0 )) have to be set. The RLSSVF estimates obtained with OSVF and ( θ(t 0 ), P osvf (t 0 )) are identical to the RLSSVF estimates obtained with NSVF and ( θ(t 0 ), P nsvf (t 0 )).

Notice that the ordinary version could introduce a significant signal scaling, leading to numerical problems, that deteriorate the identification results1 . Suggested choice/procedure 3.1 In order to avoid signal scaling problems we suggest to use in general the normalized version of the SVF.

Choice of λ svf

In [Canudas de [START_REF] De | Recursive estimation of the continuous-time process parameters[END_REF], it has been shown that the parameter convergence of RLSSVF depends on the cut-off frequency of the SVF λ svf , the sampling time and the input signal. Regarding the cut-off frequency of the SVF, the larger the value of λ svf the quicker the convergence rate. On the other hand, the smaller the sampling time the quicker the convergence rate.

The cut-off frequency λ svf is a user parameter that should be chosen somewhat larger than the system bandwidth. In the LTV case and especially for systems with relatively large variations of the bandwidth, the specification of λ svf can be critical since the system bandwidth is time-varying. For its choice, the convergence aspects discussed above should be also keep in mind.

Suggested choice/procedure 3.2 In general, the identification of LTI should be easier than the identification of LTV models. Thus, assuming that the parameters are slowly varying, one approach to have an idea about λ svf is to estimate local LTI models on segments of the training data. One option is to use SRIVC.

Digital implementation of CT filtering operations

The digital implementation issues of the CT filtering operations using an SVF (see (1.20) or (1.21)) are well-known in CT model identification. They should be treated in an appropriate way since errors generated by the digital implementation can have an impact on the quality of the estimated model.

When the intersample behavior of the filter input is piecewise constant, i.e. a zeroorder hold assumption is considered, an exact solution to the filtering operation can be obtained. This is valid both for off-line and on-line computations. When the filter input is not piecewise constant, a traditional assumption is to consider it to be linear. Then, we could use a first-order hold assumption or the bilinear (or Tustin) transform. An advantage of the latter is that it maps a stable region in the s-plane into a stable region in the z-plane [Franklin et al., 1998, p. 193]. Notice that the accuracy of these approaches depend on the sampling time, namely the faster the sampling the more accurate the filter output.

In this thesis, when the intersample behavior is assumed to be linear, the Tustin transform is considered.

Remark 3.1 In the off-line case, when the intersample behavior is assumed to be linear, a first-order hold assumption can be used without introducing any delay in the signal. However, in the on-line case, both the first-order hold assumption and the bilinear transform introduce a delay.

The digital implementation of the CT filtering operations using an SVF is carried out from a time-invariant state space representation. In Appendix A, a general digital implementation considering a time-varying state space model is presented. That holds also for the LTI SVF but with constant matrices.

Suggested choice/procedure 3.3 When the filter input is not piecewise constant, we suggest to use the bilinear transform. In order to obtain accurate filter outputs, the sampling frequency should be as high as possible.

Aspects related to the IV approaches

In this section, the discussions will be done considering only the RSRIVC method, although they will hold for RIVSVF, CLRSRIVC and CLRIVSVF as well.

Start of the IV methods

RSRIVC is initialized using RLSSVF. The switch between these algorithms takes place at some time instant t s,iv , which is when RSRIVC defined in (2.29) starts to operate. The prefiltered time-derivatives and instruments are computed from the final RLSSVF estimate; afterwards, previous RSRIVC estimates are used.

The switch to (2.29) at t s,iv should be done once convergence of the RLSSVF estimates is achieved. The convergence of RLSSVF could be automatically determined based on checks of variability or steady state condition of:

• ∆ θ(t k ) = θ(t k-1 ) -θ(t k )
• the size of the parameter covariance matrix P (t k ), which can be measured through its trace or determinant

• the condition number of the parameter covariance matrix P (t k )

• the prediction error ε(t k )

• the output error ε y (t k ) defined in (1.60)

For instance, the variability of ∆ θ(t k ) can be measured through its variance that may be computed recursively with an exponential moving average [START_REF] Bittencourt | An algorithm for finding process identification intervals from normal operating data[END_REF]. The steady state condition, for instance of the output error, can be assessed using the approach proposed in [START_REF] Cao | An efficient method for on-line identification of steady state[END_REF]. Note that the automatic detection of convergence of RLSSVF requires to set threshold values, which depend on the system.

Suggested choice/procedure 3.4 A simple option to set t s,iv is to use RLSSVF with some training data, and to check convergence based on the prediction error or the output error. For the latter, the frequency content is lower, and therefore it might be easier to visualize convergence in a plot.

Stability test

In off-line estimation using SRIVC, the prefilter (2.9) and the auxiliary model (2.6), depend on the estimates of the previous iteration. To ensure that both the prefilter and auxiliary model are stable, a stability test is required (see e.g. [START_REF] Garnier | Direct continuous-time approaches to system identification. Overview and benefits for practical applications[END_REF]). Analogously, in on-line estimation using RSRIVC, the prefilter (2.19) and the auxiliary model (2.16) depend on the estimates of the previous recursion θ(t k-1 ). Then, a stability check is also required. However, in this case it is more complicated, since instead of dealing with an LTI system, we have an LTV system. Both the prefilter and the auxiliary model can be written as an LTV state-space model

M ss ẋ(t) = F (t)x(t) + G(t)u(t) y(t) = H(t)x(t) + J(t)u(t) (3.7)
Let us recall first how to check stability of (3.7). This can be done using the following lemma (see [Rugh, 1996, p. 206])

Lemma 1 Suppose the system (3.7) is uniform exponentially stable, and there exist finite constant µ 1 , µ 2 and µ 3 such that for all t

G(t) ≤ µ 1 , H(t) ≤ µ 2 , J(t) ≤ µ 3 (3.8)
Then the state equation also is uniformly bounded-input, bounded-output stable.

More details are given in Appendix B, where we recall some important theorems about LTV system stability. Assuming that G(t), H(t) and J(t) are bounded, i.e. (3.8) is fulfilled, it remains to check uniform exponential stability. That can be done in practice by solving a linear matrix inequality problem.

A more simple approach is the following. In general, uniform exponential stability cannot be characterized by the location of the eigenvalues of the matrix F (t). However, if it is assumed that the parameters vary slowly, uniform exponential stability is guaranteed if the eigenvalues of F (t) are in the left half plane [START_REF] Ilchmann | Sufficient conditions for stability of linear time-varying systems[END_REF]. Such procedure is usually considered in recursive estimation (see e.g. [Ljung and Söderström, 1983, p. 93]) and it is also the one adopted in this thesis. Now we recall a stability test for θ(t k ) [START_REF] Ljung | System Identification. Theory for the User[END_REF]. Let us define first a set D s containing stable estimates, i.e. D s = {θ | M ss has all eigenvalues in the left half plane} (3.9)

If we denote the estimate before the stability check by θ⋆ (t k ), and after it by θ(t k ), then

θ(t k ) = θ⋆ (t k ) if θ⋆ (t k ) ∈ D s θ⋄ (t k ) if θ⋆ (t k ) / ∈ D s (3.10)
where θ⋄ (t k ) is a stable estimate obtained

• from the previous recursion, i.e. θ(t k-1 ), or

• by reflecting the eigenvalues of F (t) (obtained from θ⋆ (t k )) into the left half plane. This means that the sign of the real part of the unstable eigenvalues is changed.

Correlation between filtered instruments and filtered noise

In the batch estimation of LTI systems, consistent estimates can be obtained with SRIVC if the two conditions specified in Section 2.2.1 are fulfilled. The condition C2 is given by Ē{F

(p, θ o )ζ(t k , θ o )F (p, θ o )v o (t k , θ o )} = 0 or Ē{ζ f (t k , θ o )v o,f (t k , θ o )} = 0 which means that ζ f (t k ) and v o,f (t k ) are uncorrelated.
In the on-line case with RSRIVC, C2 takes the form (t k ). Thus, (3.11) will be slightly violated. This is analogous to what happens in DT identification when IV techniques are used (see [Ljung and Söderström, 1983, p. 237] and also [Young, 1984, p. 132]).

Ē{F (p, θ(t k-1 ))ζ(t k , θ(t k-1 ))F (p, θ(t k-1 ))v o (t k )} = Ē{ζ f (t k , θ(t k-1 ))v o,f (t k , θ(t k-1 ))} = 0 (3.11)
To reduce this correlation, as initially suggested in [Young, 1984, p. 132], the estimates can be filtered. The DT filter, denoted here by M θ (q -1 ) has the additional effect of smoothing the estimates, which improves the stability of the algorithm [Ljung and Söderström, 1983, p. 315]. A scheme of RSRIVC, including the filter M θ (q -1 ) and the stability test described previously, is shown in Figure 3.2. The filter M θ (q -1 ) can be an infinite impulse response (IIR) or a finite impulse response (FIR) filter.

Suggested choice/procedure 3.5 As proposed in [START_REF] Wong | Identification of linear discrete time systems using the instrumental variable method[END_REF], to reduce correlation between filtered instruments and filtered noise, we consider as default option for M θ (q -1 ) a delay, i.e. 2M θ (q -1 ) = q -1 An alternative to filtering the estimates is to use regularization. This is discussed in Section 3.6.2.

Estimation of the normalized covariance matrix Q n

To apply Kalman filter based methods, the user has to specify the normalized covariance matrix Q n . The performance of these methods, in terms of tracking ability and noise sensitivity, depends on Q n . The maximum likelihood method is one of the available approaches that can be used to estimate this hyperparameter Q n [START_REF] Young | Recursive Estimation and Time-series Analysis. An Introduction for the Student and Practitioner[END_REF].

Review of maximum likelihood method

The maximum likelihood method is presented in detail in [Åstrom, 1980] (see also [Söderström andStoica, 1989, Ljung, 1999]. Here we recall this method considering a simple Figure 3.2: Diagram of RSRIVC including a DT filter M θ (q -1 ) and a stability test.

case. We assume that the system is modeled by a random walk for the parameters and a linear regression, which could represent a LTV ARX model, either in CT or DT, i.e.

M θ(t

k ) = θ(t k-1 ) + w(t k ) z(t k ) = ψ T (t k )θ(t k ) + e(t k ) (3.12)
where w(t k ) and e(t k ) are independent zero-mean DT Gaussian noise processes with covariance matrix Q w variance σ 2 e , respectively. For the estimation of (3.12), there is available N measurements of the regressor and the output denoted by Ψ N = {ψ(t k )} N k=1 and Z N = {z(t k )} N k=1 , respectively. Assuming that Z N is a random sequence whose observations are not independent, the likelihood function, i.e. the probability density function of the observations Z N given by the parameters θ is defined by

p θ (Z N |Ψ N ) = p(Z N |Ψ N , θ) = N k=1 p(z(t k )|Ψ k , θ) (3.13)
The maximum likelihood method consists in finding the parameters θ such that the measurements Z N are as likely as possible [START_REF] Ljung | System Identification. Theory for the User[END_REF]. Mathematically, this is expressed as follows θ = arg max

θ p θ (Z N |Ψ N ) (3.14)
Due to numerical reasons, the logarithm of the likelihood, called log-likelihood, is usually maximized instead. Then, the maximization of the log-likelihood, denoted by l(θ), yields the maximum likelihood estimates, i.e.

θ = arg max θ log p θ (Z N |Ψ N ) = arg max θ l(θ) (3.15)
If it is assumed that e(t k ) is a zero-mean Gaussian noise sequence with variance σ 2 e , then, the conditional densities in (3.13) are Gaussian, and the likelihood is given by

p θ (Z N |Ψ N ) = N k=1 1 2πr(t k ) exp - [z(t k ) -ẑ(t k )] 2 2r(t k ) (3.16)
with conditional mean ẑ(t k ) and conditional variance r(t k ) that are computed with the Kalman filter, i.e.

ẑ(t k ) = Ē[z(t k )|Ψ k , θ] = ψ T (t k ) θ(t k |t k-1 ) (3.17) r(t k ) = var[z(t k )|Ψ k , θ] = σ 2 e + ψ T (t k )P * (t k |t k-1 )ψ(t k ) (3.18)
Notice that in (3.18) the standard form of the Kalman filter is used, and P * (t k |t k-1 )/σ 2 e = P (t k |t k-1 ). Then, the log-likelihood is given by

l(θ) = - 1 2 N log 2π + N k=1 log r(t k ) + N k=1 [z(t k ) -ẑ(t k )] 2 r(t k ) (3.19)
Finally, the maximization of (3.19) is equivalent to

θ = arg min θ N k=1 log r(t k ) + N k=1 [z(t k ) -ẑ(t k )] 2 r(t k ) (3.20)

Estimation of Q n using maximum likelihood

In this section we will show how to estimate the normalized covariance matrix in a recursive off-line fashion considering the CT LTV OE model (2.28). This model can be written as the filtered linear regression (2.24) plus the random walk model for the parameters, i.e.

M θ(t k ) = θ(t k-1 ) + w(t k ) y (na) f (t k ) = ϕ f T (t k )θ(t k ) + v f (t k ) (3.21)
Notice that (3.21) is similar to (3.12). The difference is that the filtered noise v f (t k ) in (3.21) is in general not Gaussian, although e(t k ) in (2.28) is assumed to be Gaussian. However, when RSRIVC is used, v f (t k ) will be approximately asymptotically Gaussian. Thus, considering that the estimation is performed with RSRIVC, we assume that v f (t k ) is Gaussian and the maximum likelihood method presented in the previous section can be applied to estimate the normalized covariance matrix Q n of the Kalman filter. We will consider that Q n is diagonal and its elements are denoted by q n,i with i = 1, . . . , d. Therefore, Q n is estimated by solving an optimization problem similar to (3.20) but with a parameter vector η that includes the constant scalars q n,i , i.e.

η = q n,1 . . . q n,d (3.22) Then, η = arg min η V = arg min η N k=1 log r(t k ) + N k=1 [z(t k ) -ẑ(t k )] 2 r(t k ) (3.23)
The conditional mean

ẑ(t k ) = ŷ(na) f (t k ) = ϕ f T (t k ) θ(t k |t k-1 ) (3.24)
and conditional variance

r(t k ) = σ 2 e + ϕ f T (t k )P * (t k |t k-1 )ϕ f (t k ) (3.25)
are computed with RSRIVC. The procedure requires to set initial values for Q n , that should be defined considering Remark 1.3. We should also keep in mind that the user has to specify P (t 0 ), θ(t 0 ) and λ svf . Additionally it might by necessary to incorporate bounds (l b,i , u b,i ) for q n,i , as suggested in [START_REF] Bavdekar | Identification of process and measurement noise covariance for state and parameter estimation using extended Kalman filter[END_REF] for a similar problem.

Suggested choice/procedure 3.6 From a physical insight of the system, it might be known that some parameters are constant or that some parameters vary at the same rate. To simplify the estimation of Q n through maximum likelihood, it is important to use that prior knowledge. On the other hand, since we assume that the parameters are slowly varying, it is possible to identify local LTI models on segments of some training data. Then, from the set of estimated LTI model, we can learn about the parameter variations.

Numerical example

For the following example, the normalized covariance matrix is estimated using maximum likelihood. The data generating system is given by

S p + a o 1 (t) x(t) = b o 0 u(t) y(t k ) = x(t k ) + e(t k ) (3.26)
with b o 0 = 3 and a o 1 (t) = 2 + 0.8t/250 The sampling time is set to 0.025 s and the total simulation time is 250 s. The input is a PRBS and the DT measurement noise is a zero-mean, Gaussian noise with constant variance 3.0 • 10 -3 . Part of the input-output data is shown in Figure 3.3.

It is assumed that the parameter b 0 is known to be constant, and thus q n,2 = 0. Therefore, η is the scalar q n,1 corresponding to a 1 . The optimization problem (3.23), with initial value q n,1 = 10 -5 , lower bound l b = 10 -6 , upper bound u b = 10 -4 and λ svf = 3 (rad/s), yields q n,1 = 3.3 • 10 -5 . The plot cost function V vs η is shown in Figure 3.4.

The normalized covariance matrix is then given by

Q n = diag{3.3 • 10 -5 0} (3.27)
The true parameters are compared with the RSRIVC estimates computed with (3.27) in Figure 3.5.

Handling of input/output initial conditions

In the identification of DT models, the predictor depends on unknown input/output initial conditions. To circumvent that problem, the initial conditions can be assumed to be zero [Ljung, 1999, p. 320]. Analogously, in the derivation of the LSSVF method in Section 1.2.2, the initial conditions of the LTI system (1.1) have been ignored. That is, (1.22) has been obtained assuming that the initial conditions are zero. These initial conditions can be estimated together with the model parameters, as it was shown in [START_REF] Garnier | A time-domain approach to continuous-time model identification of highly resonant wide-band systems[END_REF]. The goal of this section is to illustrate the impact of the initial conditions in the recursive estimates. To this end, we derive the RLSSVF for the LTI case considering initial conditions.

Problem formulation

We start by formulating the problem again. Let us simplify the problem by assuming that there is no measurement noise in the data generating system, which is given by

S : A o (p)y(t) = B o (p)u(t) (3.28)
with polynomials A o (p) and B o (p) relatively coprime. The system is assumed to be subject to an arbitrary set of input/output initial conditions

u 0 = u(0) u (1) (0) • • • u (n b -1) (0) , (3.29) y 0 = y(0) y (1) (0) • • • y (na-1) (0) . (3.30)
Moreover, it is assumed that the system belongs to the model set defined by

M : A(p, θ)y(t) = B(p, θ)u(t) (3.31)
A(p, θ) and B(p, θ) are the following polynomials

B(p, θ) = b 0 p n b + b 1 p n b -1 + . . . + b n b A(p, θ) = p na + a 1 p na-1 + . . . + a na
where n a ≥ n b . Then, the identification problem is to recursively estimate the parameters that characterize the model structure given by (3.31), based on the sequences {u(t k ), y(t k )} N ′ k=1 , where N ′ is the number of samples which increases by one with every recursion.

Recursive least squares state-variable filter method

Consider the Laplace transform of the differential equation defined in (3.31),

s na Y (s) + na i=1 a i s na-i Y (s) = n b i=0 b i s n b -i U (s) + na-1 i=0 c i s na-i (3.32) or A(s)Y (s) = B(s)U (s) + C(s) (3.33)
where s represents the Laplace variable while Y (s) and U (s) are the Laplace transforms of y(t) and u(t) respectively. The coefficients c i depend on the unknown parameters a i and b i as well as the unknown initial conditions. Notice that if there was measurement noise, the coefficients c i would depend on initial conditions of the noise as well.

Applying a filter F (s) to (3.33) yields

A(s)F (s)Y (s) = B(s)F (s)U (s) + F (s)C(s), (3.34)
The SVF has the following form

F (s) = λ svf s + λ svf na (3.35)
where λ svf is the cut-off frequency of the filter.

Let us define a filter F i (s) to be

F i (s) = (λ svf ) na s i (s + λ svf ) na for i = 1, . . . , n a (3.36)
Then, (3.34) in terms of F i (s) takes the following form

F na (s) + a 1 F na-1 (s) + . . . + a na F 0 (s) Y (s) = b 0 F n b (s) + . . . + b n b F 0 (s) U (s) + c 1 F na (s) + . . . + c na F 0 (s) (3.37)
In terms of time-domain signals, (3.37) can be written as

[F na y](t) + a 1 [F na-1 y](t) + . . . + a na [F 0 y](t) = b 0 [F n b u](t) + . . . + b n b [F 0 u](t) + c 1 f na-1 (t) + . . . + c na f 0 (t) (3.38)
where

[F i y](t) = y (i) f (t) = f i (t) * y(t) [F i u](t) = u (i) f (t) = f i (t) * u(t)
with f i the corresponding impulse response of (3.35) and * denoting the convolution operator. The filter outputs [F i y] and [F i u] are prefiltered time-derivatives of the outputs and inputs in the bandwidth of interest. Then, (3.38) can be expressed as the linear regression

y (na) f (t k ) = ϕ f T (t k )θ (3.39)
The parameter vector is given by

θ = a 1 . . . a na b 0 . . . b n b c 1 . . . c na T (3.40)
and the regressor by

ϕ f (t k ) =     ϕ f,y (t k ) ϕ f,u (t k ) ϕ ic (t k )     (3.41) with ϕ T f,y (t k ) = -y (na-1) f (t k ) y (na-2) f (t k ) • • • y f (t k ) ϕ T f,u (t k ) = u (n b ) f (t k ) u (n b -1) f (t k ) • • • u f (t k ) ϕ T ic (t k ) = f na-1 (t k ) f na-2 (t k ) • • • f 0 (t k )
From (3.39) it is straightforward to estimate θ using the RLSSVF-FF algorithm (1.57) with λ = 1 (for the LTI case considered here). The estimation of initial conditions implies an extra computational burden. On the other hand, neglecting its effect can lead to overshoots in the estimates. In order to avoid these overshoots, another option is to start the estimation after transients due to initial conditions have died out. Suggested choice/procedure 3.7 The proposed identification strategy with RLSSVF is the following: prefiltered time-derivatives are computed from t = 0, but estimates are computed only from the time instant when the effect of the initial conditions are negligible. This time instant when the least squares approach starts to operate is denoted by t s,ls . In the following example, these ideas are illustrated and further discussed.

Remark 3.2 An alternative to deal with overshoots is to use regularization based methods, as pointed out e.g. in [START_REF] Stoica | Exact initialization of the recursive least squares[END_REF]. The use of regularization is discussed in Section 3.6.

Numerical example

The data generating system is noise-free and given by

(p + a o 1 )y(t) = b o 0 u(t) (3.42)
with a o 1 = 3 and b o 0 = 30. The sampling time is T s = 0.01 s and the input is a square wave. In the input-output data shown in Figure 3.6, the transient response has been removed. Therefore, the output is periodic for t ≥ 0 and the period is 2 s. The initial condition is y(0) = 18.1.

We are interested in identifying a CT model that represents (3.42) considering initial conditions. For this case, we have a linear regression as in (3.39) with

ϕ f T (t k ) = -y f (t k ) u f (t k ) f 0 (t k ) θ T = a 1 b 0 c 1 where f 0 (t k ) is the impulse response corresponding to the SVF F 0 (s) = λ svf s + λ svf (3.43)
Notice that c 1 is here an estimate of y(0). The parameters are estimated with RLSSVF using λ svf = 4 rad/s. In Figure 3.7, we compare RLSSVF estimates obtained with and without taking into account the initial conditions. We can see that by taking into account the initial conditions the results are better in terms of convergence rate and bias. The estimate ĉ1 (t k ) is very close to the true value (y(0) = 18.1), as it can be seen in Figure 3.8, where we plot additionally the impulse response f 0 (t k ). Now let us take a look at the overshoots in the estimates that cannot be seen in Figure 3.7 because of the chosen axes. The overshoots are defined as follows

overshoot = max[ θ(t k )] -θ o θ o × 100 (3.44)
The results, presented in Table 3.1, show that large overshoots are obtained when the initial conditions are not considered. By including the estimate of c 1 , the overshoots are significantly reduced.

As it was mentioned in the previous section, another option is to start the estimation after transients due to initial conditions have died out. From Figure 3.8, for the chosen value of λ svf , we can see that the effect of the initial condition is negligible for t ≥ 1 s. The RLSSVF estimates with t s,ls = 2 s is shown in Figure 3.9. Note that since the input-output data is periodic and there is no noise, the observations for t = 2 s are the same than for t = 0. These simulation conditions have been deliberately chosen in order to show only the effect of the initial conditions. In Table 3.1, we can see that the overshoots are the smallest in this case. From (3.43) it can be seen that the effect of the initial conditions is weighted by λ svf . Therefore, if RLSSVF estimates are computed without initial conditions, the larger the value of λ svf , the larger the overshoots will be. This has been also pointed out in [START_REF] Johansson | Identification of continuous-time models[END_REF], for an approach which is similar to RLSSVF.

Numerical issues

In off-line identification, least squares estimates can be computed using the normal equations. Nonetheless, the conditioning can be improved by using QR factorization. Another complementary method to improve conditioning is signal scaling. While off-line estimation is usually done in a desktop computer, on-line estimation is meant to be performed in a microcontroller, where numerical aspects are more important. Roundoff errors introduced by the microcontroller due to its poor machine precision might yield a negative definite covariance matrix. In recursive estimation, QR factorization can be also implemented in a recursive fashion with forgetting factor (see [START_REF] Bittencourt | An algorithm for finding process identification intervals from normal operating data[END_REF]). Thus, the forgetting factor based RLSSVF method (1.57) can be implemented using QR factorization. However, the Kalman filter based RLSSVF cannot be handled in this way. A solution is Potter's algorithm which corresponds to a square root filtering method (see e.g. [START_REF] Niedźwiecki | Identification of Time-varying Processes[END_REF][START_REF] Grewal | Kalman Filtering[END_REF][START_REF] Simon | Optimal State Estimation[END_REF]). For Potter's algorithm, there are both a Kalman filter based approach and a forgetting factor based approach.

Next, we discuss first the numerical issues in off-line identification. After showing how to measure conditioning in LSSVF, we present the signal scaling approach for LSSVF to improve conditioning. Then, the numerical issues in recursive identification are addressed. First we discuss solutions for LS based methods and at the end for the IV based methods developed in Chapter 2. In Section 4.1, the proposed approaches to improve conditioning are tested with real data.

Conditioning in LSSVF

For LSSVF presented in Section 1.2.2, let us consider the minimization problem (1.10) in matrix form. Then, the cost function is given by

V = 1 N Y -Φθ 2 (3.45)
where

Y =        y (na) f (t 1 ) y (na) f (t 2 ) . . . y (na) f (t N )        , Φ = Φ y Φ u with Φ y = -         y (na-1) f (t 1 ) y (na-2) f (t 1 ) • • • y f (t 1 ) y (na-1) f (t 2 ) y (na-2) f (t 2 ) . . . . . . . . . . . . y (na-1) f (t N ) y (na-2) f (t N ) • • • y f (t N )         Φ u =         u (n b ) f (t 1 ) u (n b -1) f (t 1 ) • • • u f (t 1 ) u (n b ) f (t 2 ) u (n b -1) f (t 2 ) . . . . . . . . . . . . u (n b ) f (t N ) u (n b -1) f (t N ) • • • u f (t N )        
The optimization problem with objective function (1.10) corresponds to the normal equations

Φ T Φ θ = Φ T Y (3.46) or R N θ = f N (3.47) which have a closed-form solution θ = R -1 N • f N (3.48)
There are different ways to check whether R N in (3.48) is well-conditioned. If R N is well-conditioned, it is invertible. The invertibility of a matrix can be evaluated by computing its determinant which must be different than zero. However, this might not be sufficient to detect numerically ill-conditioned matrices whose determinant might be close to zero [START_REF] Shardt | Data quality assessment of routine operating data for process identification[END_REF]]. An alternative is to use the condition number, which is defined by [Quarteroni et al., 2000, p. 60]

κ p (R N ) = R N p R -1 N p (3.49)
where • p is the p-norm. For the infinity norm (p = ∞), the values tend to be too conservative, and usually 2-norm (p = 2) is preferred [START_REF] Shardt | Data quality assessment of routine operating data for process identification[END_REF]. The 2-norm condition number is defined by

κ 2 = σ max σ min (3.50)
where σ max , σ min is the maximum and minimum singular values, respectively. For a symmetric positive definite matrix (like in this case for R N )

κ 2 = λ max λ min (3.51)
with λ max , λ min the maximum and minimum eigenvalues, respectively [Quarteroni et al., 2000, p. 61].

Remark 3.3

The condition number κ 2 (M ) lies in the interval [1 ∞) and the optimal conditioning is 1, which is obtained in DT identification with a FIR model with white noise as input [START_REF] Van Den Hof | System identification with generalized orthonormal basis functions[END_REF].

Signal scaling in LSSVF

Signal scaling can be used to improve conditioning of the normal equations. In the offline identification of a DT linear regression model using least squares, the regressor is defined in term of past outputs and inputs (see (1.5)). Thus, the linear regression can be readily scaled by changing the magnitude of the input-output signals. The situation is more tricky in the off-line identification of a CT linear regression model using LSSVF. In that case, we arrive at a filtered linear regression (see (1.24)) whose regressor is defined in terms of prefiltered time-derivatives of the input-output signals. Thus, the filtered linear regression cannot be longer scaled by changing the magnitude of the input-output signals.

Next we show how to apply signal scaling in LSSVF. The cost function (3.45) can be modified as follows

V = 1 N Y -ΦM s M -1 s θ 2 (3.52)
where M s is a diagonal matrix that is used to obtain a scaled solution. The matrix M s is defined by

M s =     m s,1 0 0 0 . . . 0 0 0 m s,n θ     (3.53) with m s,i = 1 Φ i (3.54)
and Φ i the column vectors that constitute Φ, i.e.

Φ = Φ 1 . . . Φ n θ (3.55) Let us define θs = M -1 s θ (3.56a) Φ s = ΦM s (3.56b)
Then, the cost function (3.52) may by written as

V = 1 N Y -Φ s θ s 2 (3.57)
Thus, the scaled LSSVF estimate that minimizes (3.57) is given by,

θs = [Φ T s Φ s ] -1 Φ T s Y (3.58)
In general, the filtered noise v f (t k ) is not Gaussian. However, if the polynomial A(p, θ) in 1.2 and the polynomial E(p) in (1.19) are equal, v f (t k ) is Gaussian and the parameter covariance matrix can be estimated by

P = σ2 v f [Φ T Φ] -1 (3.59)
with σv f (t k ) the sample variance of v f (t k ). Equation (3.59), that will allow us to introduce the signal scaling approach in the recursive case in the next section, can be written as follows

P [Φ T Φ] = σ2 v f (3.60) Replacing (3.56b) in (3.60) yields P M -1 s Φ T s Φ s M -1 s = σ2 v f (3.61)
Multiplying (3.61) from the left by M -1 s we obtain

M -1 s P M -1 s Φ T s Φ s M -1 s = σ2 v f M -1 s (3.62)
Let us define

P s = M -1 s P M -1 s (3.63)
Then, multiplying (3.62) from the right by M s yields

P s = σ2 v f [Φ T s Φ s ] -1 (3.64)
Then, by choosing M s according to (3.53), κ 2 (P s ) < κ 2 (P ). From the properties of the condition number, note that for R N defined in (3.48), κ 2 (R N ) = κ 2 (P ).

Signal scaling in RLSSVF

Signal scaling can be also applied in recursive estimation. The expression (3.56a) shows us how to scale θ(t k ) in the KF based RLSSVF algorithm defined in (1.59). Note that although (3.63) is in general not valid, it also shows us how to modify P (t k ) in (1.59). Then, in analogy to the off-line case (see (3.56a) and (3.63)), for the time instant t k , we can write

θs (t k ) = M -1 s θ(t k ) (3.65a) P s (t k ) = M -1 s P (t k )M -1 s (3.65b)
Replacing (3.65) in (1.59) we obtain the scaled version of the KF based RLSSVF:

Prediction step:

θs (t k |t k-1 ) = θs (t k-1 ) (3.66a) P s (t k |t k-1 ) = P s (t k-1 ) + Q n,s (3.66b)
Correction step:

θs (t k ) = θs (t k |t k-1 ) + M -1 s L(t k )ε(t k ) (3.66c) ε(t k ) = y (na) f (t k ) -ϕ f T (t k )M s θs (t k |t k-1 ) (3.66d) L(t k ) = M s P s (t k |t k-1 )M s ϕ f (t k ) 1 + ϕ f T (t k )M s P s (t k |t k-1 )M s ϕ f (t k ) (3.66e) P s (t k ) = M -1 s [I -L(t k )ϕ f T (t k )]M s P s (t k |t k-1 ) (3.66f)
where

Q n,s = M -1 s Q n M -1 s .
Notice that in this recursive version M s could be computed with some representative training data. It is clear from (3.66) that the disadvantage of this approach are the additional operations that are introduced.

Review of Potter's algorithm

Next we recall Potter's algorithms which have been developed for DT estimation. Since the only difference in the CT case w.r.t. the DT case is how the linear regression is built, we present directly the CT version of the approaches.

The conditioning of the recursive estimation can be also improved by using square root filtering algorithms (see e.g. [START_REF] Grewal | Kalman Filtering[END_REF]). Instead of propagating the covariance matrix P (t k ), they propagate S(t k ) which is the square root of P (t k ), i.e.

P (t k ) = S(t k )S T (t k ) (3.67)
The benefit of this approach is clear, since the condition number of S(t k ) is the square root of the condition number of P (t k ), i.e.

κ 2 (P (t k )) = κ 2 2 (S(t k )) (3.68)
One of these approaches is Potter's algorithm, which consists on writing the covariance matrix of the correction step as follows

S(t k )S T (t k ) = S(t k-1 ) I d - f (t k )f T (t k ) β(t k ) S(t k-1 ) (3.69)
where

f (t k ) = S T (t k-1 )ϕ f (t k ) (3.70a) β(t k ) = 1 + f T (t k )f (t k ) (3.70b)
It has been shown that the square root of P (t k ) can be updated as follows (see e.g. [START_REF] Niedźwiecki | Identification of Time-varying Processes[END_REF])

S(t k ) = S(t k-1 ) I d - f (t k )f T (t k ) α(t k ) (3.71)
where

α(t k ) = β(t k ) + β(t k ) (3.72)
The Kalman filter based Potter's algorithm is given by Prediction step:

θ(t k |t k-1 ) = θ(t k-1 ) (3.73a) S T (t k |t k-1 ) 0 = T S T (t k-1 ) Q T /2 n (3.73b)
Correction step:

θ(t k ) = θ(t k |t k-1 ) + L(t k )ε(t k ) (3.73c) ε(t k ) = y (na) f (t k ) -ϕ T f (t k ) θ(t k |t k-1 ) (3.73d) f (t k ) = S T (t k |t k-1 )ϕ f (t k ) (3.73e) β(t k ) = 1 + f T (t k )f (t k ) (3.73f) α(t k ) = β(t k ) + β(t k ) (3.73g) L(t k ) = S(t k |t k-1 )f (t k ) β(t k ) (3.73h) S(t k ) = S(t k |t k-1 ) I d - f (t k )f T (t k ) α(t k ) (3.73i)
where T is a 2n θ × 2n θ orthogonal matrix that can be found using Givens rotations, Householder transformation, etc. (see more in [Simon, 2006, p. 162]). Regarding the matrix S(t k |t k-1 ), note that it is not unique. Algorithm (3.73) will be called PotterSVF.

Signal scaling can also be introduced in PotterSVF as it is shown in Appendix C.

Signal scaling in RSRIVC

The RSRIVC method presented in Section 2.3.3 is the on-line counterpart of SRIVC.

The SRIVC estimate can be written as in (3.48), i.e.

θi = N k=1 ζ f (t k , θi-1 )ϕ f T (t k , θi-1 ) -1 • N k=1 ζ f (t k , θi-1 )y (na) f (t k , θi-1 )
This approach corresponds to the so-called non-symmetric version of the IV technique [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], since in this case, the matrix in the first square brackets is non-symmetric. An alternative is the symmetric SRIVC defined by θi =

N k=1 ζ f (t k , θi-1 )ζ T f (t k , θi-1 ) -1 • N k=1 ζ f (t k , θi-1 )y (na) f (t k , θi-1 )
where in this case, the matrix in the first square brackets is symmetric. Analogously, we can derive the recursive counterpart of the symmetric SRIVC.

The drawback of the square root filtering method is that it cannot be applied to nonsymmetric IV techniques, since the corresponding covariance matrix is non-symmetric and it cannot be written as (3.67). Such an issue could be circumvented by using the symmetric IV approach. However, it is known that the non-symmetric IV approach gives better results than the symmetric IV approach [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]. A solution to improve conditioning in the proposed recursive IV methods presented in Chapter 2 is then to use signal scaling. Analogously to RLSSVF (see Section 3.5.3), the scaled version of the KF based RSRIVC method is given by Prediction step:

θs (t k |t k-1 ) = θs (t k-1 ) (3.74a) P s (t k |t k-1 ) = P s (t k-1 ) + Q n,s (3.74b)
Correction step:

θs (t k ) = θs (t k |t k-1 ) + M -1 s L(t k )ε(t k ) (3.74c) ε(t k ) = y (na) f (t k ) -ϕ f T (t k )M s θs (t k |t k-1 ) (3.74d) L(t k ) = M s P s (t k |t k-1 )M s ζ f (t k ) 1 + ϕ f T (t k )M s P s (t k |t k-1 )M s ζ f (t k ) (3.74e) P s (t k ) = M -1 s [I -L(t k )ϕ f T (t k )]M s P s (t k |t k-1 ) (3.74f)
where

Q n,s = M -1 s Q n M -1 s .

Anti-windup techniques

In system identification, in order to obtain accurate models, the input u of the system has to be persistently exciting (or sufficiently rich). Loosely speaking this means that the input should excite all modes of the system [START_REF] Ljung | Theory and Practice of Recursive Identification[END_REF]. In formal terms, for off-line estimation using LSSVF, the input is persistently exciting if the matrix R N in (3.48) is non-singular. Similarly, in on-line estimation using RLSSVF, the persistence of excitation condition is defined by

i+s k=i ϕ f (t k )ϕ f T (t k ) c • I ∀k (3.75)
with s and c positive constants. Notice that for stochastic regressors, the deterministic persistence of excitation condition (3.75) is very difficult to fulfill (see [START_REF] Ljung | Adaptation and tracking in system identification -a survey[END_REF]] and [Niedźwiecki, 2000, p. 72]). In our case, if the input is deterministic like a PRBS, and the output is contaminated with stochastic noise, the regressor will have both deterministic and stochastic components. Then (3.75) could be used to check the persistency of excitation. For the recursive IV based methods, the counterpart of (3.75) is given by

i+s k=i ζ f (t k )ϕ f T (t k ) c • I ∀k (3.76)

Anti-windup techniques in discrete-time identification

In DT recursive identification, when the input is not sufficiently rich, recursive estimation approaches could suffer from a phenomenon called estimator windup, where the covariance matrix blow up leading to wrong estimates. To illustrate the problem, let us consider an extreme case where we want to estimate a FIR model,

y(t k ) = φ T (t k )ρ + e(t k )
with φ(t k ) = 0, ∀ t k ≥ 0. Then, the RLS-FF algorithm becomes [Åstrom and Wittenmark, 2008, p. 473] ρ(t k ) = ρ(t k-1 ) (3.77a)

P (t k ) = 1 λ P (t k-1 ) (3.77b)
The equations (3.77) are unstable and P (t k ) will grow exponentially if λ < 1. On the other hand, the KF algorithm with φ(t k ) = 0 becomes [Niedźwiecki, 2000, p. 284], [START_REF] Evestedt | Parameter and State Estimation with Information-rich Signals[END_REF] ρ(t k ) = ρ(t k-1 ) (3.78a)

P (t k ) = P (t k-1 ) + Q w (t k ) (3.78b)
The equations (3.78) are also unstable but in this case P (t k ) grows linearly.

Remark 3.4 In general, it can be expected that a linear growth is slower than an exponential growth, which means that the Kalman filter with a constant Q w could be more robust to excitation problems than the forgetting factor based approach [START_REF] Cao | Analysis of the Kalman filter based estimation algorithm: an orthogonal decomposition approach[END_REF].

Many anti-windup techniques have been proposed in the literature (see e.g. [START_REF] Salgado | Modified least squares algorithm incorporating exponential resetting and forgetting[END_REF], Parkum et al., 1992[START_REF] Gunnarsson | Combining tracking and regularization in recursive least squares identification[END_REF][START_REF] Cao | A directional forgetting algorithm based on the decomposition of the information matrix[END_REF][START_REF] Evestedt | Parameter and State Estimation with Information-rich Signals[END_REF][START_REF] Evestedt | Windup properties of recursive parameter estimation algorithms in acoustic echo cancellation[END_REF]). We can distinguish for instance the following classes:

• Directional forgetting approaches. The principle is that old data must be discarded only in the direction where new information is available. One algorithm of this class is the selective forgetting which was developed in [START_REF] Parkum | Recursive Identification of time-varying systems[END_REF], Parkum et al., 1992].

• Directional tracking approaches. The principle is to restrict tracking directions of an algorithm to the excited subspace (see [START_REF] Cao | Analysis of the Kalman filter based estimation algorithm: an orthogonal decomposition approach[END_REF]). The algorithm proposed in [START_REF] Stenlund | Avoiding windup in recursive parameter estimation[END_REF] corresponds to this class. It is a KF based algorithm and it is suitable when the parameters vary at different speeds.

• Conditional updating [Åstrom and Wittenmark, 2008, p. 477] (to be explained in Section 3.6.2).

• Regularization based approaches [START_REF] Gunnarsson | Combining tracking and regularization in recursive least squares identification[END_REF] (to be explained in Section 3.6.3).

These methods have been developed for DT estimation. Here we extend the conditional updating and regularization based approaches to the CT case using the linear filter methods that we have proposed. Since the only difference in the CT case w.r.t. the DT case is how the linear regression is built, we present directly the CT version of the approaches.

Conditional updating

This method consists in updating the estimates and covariance only when there is excitation [Åstrom and Wittenmark, 2008, p. 477]. A method that works with conditional updating is for instance the Dasgupta-Huang optimal bounding ellipsoid algorithm [START_REF] Dasgupta | Asymptotically convergent modified recursive leastsquares with data-dependent updating and forgetting factor for systems with bounded noise[END_REF]. Note that conditional updating can be done using the same techniques proposed in data quality assessment [START_REF] Peretzki | Data mining of historic data for process identification[END_REF][START_REF] Shardt | Data quality assessment of routine operating data for process identification[END_REF][START_REF] Bittencourt | An algorithm for finding process identification intervals from normal operating data[END_REF], which consists on finding suitable data for model estimation given a large dataset.

There are different measures that can be used for conditional updating and the criteria to choose one or another should be based e.g. on computational cost, reliability and sensitivity of the indexes. By sensitivity we mean how much does the measure change between rich excitation and poor excitation.

The conditional updating can be based on:

1. the conditions of persistent excitation (3.75) and (3.76) for LS and IV based methods, respectively 2. measures of the input u(t k ) or output y(t k )

3. measures of quantities coming from the recursive algorithm like the scalars ε(t k ) and ϕ f T (t k )P (t k )ϕ f (t k ) [START_REF] Åstrom | Adaptive Control[END_REF], or the matrix P (t k )

Regarding the second category, the updating could be done for instance based on a variability test of the output. In the context of data quality assessment, this is proposed in [START_REF] Bittencourt | An algorithm for finding process identification intervals from normal operating data[END_REF] by means of a recursive computation of the output variance.

Regarding the third category, the updating can be done by defining a scalar measure of the covariance matrix P (t k ). The inverse of P (t k ) can be proportional to the information matrix. Then good estimates imply a small matrix P (t k ) or a large information matrix. The size of a matrix can be measured in different ways, as it is discussed in the field of optimal experiment design [START_REF] Rojas | Robust Experiment Design[END_REF]. Possible solutions are the determinant or the trace. The other option is to compute the condition number3 of P (t k ), as it has been proposed in the studies about data quality assessment (see e.g. [START_REF] Bittencourt | An algorithm for finding process identification intervals from normal operating data[END_REF]). Since the condition number lies in the interval [1, ∞) it might be difficult to set a threshold for it. Instead, its inverse, i.e. [κ 2 (P )] -1 that lies in the interval (0, 1], can be monitored [START_REF] Bittencourt | An algorithm for finding process identification intervals from normal operating data[END_REF].

Note that in practice it is likely that conditional updating is based on several measures, in order to increase robustness of the approach. In [START_REF] Bittencourt | An algorithm for finding process identification intervals from normal operating data[END_REF] for instance, various indexes are considered for data quality assessment.

Proposed conditional updating for IV based methods

In this work, we suggest to use the recursive IV methods proposed in Chapter 2 with conditional updating based on the condition number of P (t k ). However, conditional updating based on measures provided by the algorithm, like the condition number of P (t k ), should not be done using directly an IV method. The reason is that in this case, P (t k ) depends not only on the data but also indirectly on the estimated parameters θ(t k ) through the filtered instrument ζ f (t k ). Under poor excitation, θ(t k ) is not reliable, and neither is P (t k ) as a measure for conditional updating. The solution that we propose is to run in parallel RLSSVF for the purpose of obtaining κ 2 (P ).

For RLSSVF there are two possible adaptation mechanism, the forgetting factor and the KF. The former is a better option for conditional updating because it is less robust to poor excitation (see Remark 3.4), which means that the measure κ 2 (P ) is more sensitive. Then, the conditional updating for the IV methods is based on the condition

κ 2 (P ) < γ (3.79)
where P is obtained from RLSSVF-FF and γ is the threshold which is chosen depending on the given identification problem.

When the excitation is detected as being poor such that (3.79) does not hold, the IV based algorithms must be reset to the original values, i.e. the parameter estimates and covariance matrix are set as follows θ(t k ) = θ(t 0 )

P (t k ) = P (t 0 )
Once the excitation is again rich enough, the IV approaches are initialized as usual with RLSSVF.

Levenberg-Marquardt regularization method

Regularization is a well-known technique in system identification (see e.g. [START_REF] Sjöberg | On the use of regularization in system identification[END_REF][START_REF] Pillonetto | Kernel methods in system identification, machine learning and function estimation: a survey[END_REF]). It is often used to deal with an ill-conditioned problem that can arise due to poor excitation. Anti-windup algorithms based on regularization are presented in [START_REF] Gunnarsson | Combining tracking and regularization in recursive least squares identification[END_REF][START_REF] Ch | Generalized damped least squares algorithm[END_REF][START_REF] Van Waterschoot | Optimal regularized adaptive filtering algorithms for room acoustic signal enhancement[END_REF]. It is important to keep in mind that regularization involves a bias-variance trade-off.

As it is shown in [START_REF] Sayed | A state-space approach to adaptive RLS filtering[END_REF], the forgetting factor based RLSSVF algorithm (1.57), initialized with θ 0 = θ * , P 0 = Γ -1 , corresponds to the following regularized cost function

V = N ′ k=1 λ N ′ -k [y (na) f (t k ) -ϕ f T (t k )θ] 2 + [θ -θ * ](λ (N ′ +1) Γ)[θ -θ * ] (3.80)
In (3.80) we can note that as the number of samples increases, the regularization term vanishes due to the forgetting factor. By using regularization as in (3.80) we are forcing our first estimates to be close to θ * .

Remark 3.5 To increase the convergence rates of the estimates in the transient mode, it is sometimes suggested to decrease the value of the forgetting factor (see [Ljung and Söderström, 1983, p. 279]). In this way, the tracking capability of the algorithm is increased, but at the same time it is more sensitive to noise. On the other hand, it can be observed from (3.80), that by decreasing the forgetting factor, we penalize less the regularization term, which is damping the estimates.

If we want the regularization term to remain in time, then the cost function should be defined as follows

V = N ′ k=1 λ N ′ -k [y (na) f (t k ) -ϕ f T (t k )θ] 2 + [θ -θ * ] T Γ[θ -θ * ] (3.81)
The recursive estimate that minimize (3.81) for θ * = θ(t N ′ -1 ) is presented in [START_REF] Van Waterschoot | Optimal regularized adaptive filtering algorithms for room acoustic signal enhancement[END_REF] and it is defined as follows

θ(t k ) = θ(t k-1 ) + R -1 (t k )ϕ f (t k )ε(t k ) (3.82a) R(t k ) = λR(t k-1 ) + ϕ f (t k )ϕ f (t k ) + (1 -λ)Γ (3.82b) ε(t k ) = y (na) f (t k ) -ϕ f T (t k ) θ(t k-1 ) (3.82c)
For the choice θ * = θ(t N ′ -1 ), algorithm (3.82) corresponds to the Levenberg-Marquardt regularization (LMR) approach in [START_REF] Ljung | Theory and Practice of Recursive Identification[END_REF]]. The algorithm (3.82) is written in terms of the CT linear regression that it is built with SVF. Therefore, we call it LMR-RLSSVF method, and it has two hyperparameters, λ and Γ. Notice that with this choice for θ * , we are damping the estimate such that it remains close to the previous value. The downside of LMR is that it involves matrix inversion, as it is the case for all the methods based on regularization. However, different methods have been developed to overcome or alleviate that issue (see [START_REF] Gunnarsson | Combining tracking and regularization in recursive least squares identification[END_REF] and [START_REF] Van Waterschoot | Optimal regularized adaptive filtering algorithms for room acoustic signal enhancement[END_REF]).

LMR approach and IV methods under rich excitation

While regularization introduces bias in the estimates, the aim of IV methods is to reduce it. Therefore, at first sight, it would seem that it does not make sense to use these approaches together. However, in recursive estimation of slowly time-varying models under rich excitation conditions, if we consider regularization as in (3.81) with θ * = θ(t N ′ -1 ), the bias introduced is minimum. On the other hand, in the IV techniques presented in Chapter 2, the instruments are built using previous estimates, which should vary smoothly to guarantee stability of the algorithm. A solution to get smooth estimates is to use a filter, as discussed in Section 3.2.3. Another option is to use the IV methods with regularization, i.e. the IV counterpart of algorithm (3.82). However, it is important to stress that such an approach is not suitable to deal with poor excitation, because in the considered IV methods the instruments depend on previous estimates θ(t k-1 ) that are not reliable.

Numerical example

To illustrate the windup problem and the proposed solutions, the numerical example given in Section 2.3.5 is considered under poor excitation conditions. For 0 ≤ t < 6000 s the input is a PRBS as before, and for 6000 ≤ t ≤ 9000 s the input remains equal to 1, leading to an estimator windup problem. Part of the input-output data is shown in Figure 3.10.

We present here the performance of the algorithms shown in Table 3.2. For the simulations, we use as before an SVF with cut-off frequency λ svf = 3.2 rad/s, which is slightly higher than the maximum system bandwidth. The rest of the hyperparameters has been obtained by trial and error and their values are given in Table 3.2.

Firstly, we present the results using RLSSVF-FF and RLSSVF-KF. With the former approach, good results can be expected since the parameters are varying with a similar rate. The estimates obtained using RLSSVF-FF and RLSSVF-KF are shown in Figure Figure 3.12: Condition number of P (t k ) for RLSSVF-FF and RLSSVF-KF when poor excitation appear from t ≥ 6000 s.

3.11. Under rich excitation conditions, the RLSSVF-FF estimates are similar to the ones provided by RLSSVF-KF. However this is not the case under poor excitation conditions. The condition numbers of P are plotted in Figure 3.12. We can see a faster growth of κ 2 (P ) for RLSSVF-FF; as we already mentioned, this is expected since in forgetting factor based algorithms P (t k ) may grow exponentially under poor excitation, while in KF based algorithms the growth is linear.

The estimation with RSRIVC is performed considering a conditional updating based on κ 2 (P ) of RLSSVF-FF with a threshold γ = 300 (see (3.79)). In Figure 3.13, the RSRVIC estimates are plotted. For t ≈ 6000 s, the condition κ 2 (P ) < 300 does not hold anymore, and therefore θ(t k ) = 0.

Regarding LMR-RLSSVF, the estimates are compared with the true values in Figure 3.14. Notice that under persistent excitation conditions, the RLSSVF-FF and LMR-RLSSVF estimates are similar. However, under poor excitation, κ 2 (P ) is smaller for LMR-RLSSVF, as it can be seen in Figure 3.15. Note though that κ 2 (P ) for RLSSVF-KF is even smaller (see Figure 3.12).

Conclusions

For the algorithms presented in Chapter 2, several practical aspects and implementation issues, summarized in the introduction of the chapter, have been analyzed. Next we list them together the proposed solutions. -λ svf . This is an hyperparameter of RLSSVF which is used to initialized the IV based methods. The cut-off frequency λ svf should be chosen close to the system bandwidth. One of the assumptions that we have made is that the parameters are slowly varying. Then, to have an idea about values for λ svf , we propose to estimate local LTI models on portions of some training data. Nonetheless, it is important to know that this parameter has also an impact in the parameter convergence, namely, the larger the value of λ svf , the quicker the convergence. On the other hand, in the transient mode, the larger λ svf , the larger the overshoots in the estimates can be. Through a numerical example we have shown that one way to reduce these overshoots is to start the estimation after the transients due to initial conditions have died out.

Regarding the implementation of the SVF, the ordinary version could introduce an undesired signal scaling that leads to numerical issues. Therefore we suggest to use in general the normalized version.

-t s,iv . The choice of the time instant at which RLSSVF converges and the IV methods starts to operate can be done automatically. To this end, some indexes have been presented. Based on RLSSVF estimations in some training data, we suggest to choose t s,iv according to prediction error ε(t k ) or output error ε y (t k ).

-Q n . The normalized covariance matrix of the Kalman filter can be estimated through maximum likelihood using some training data. It might be the case that it is known in advance that some parameters do not vary in time. Then, it is important that such prior information is used in the estimation through maximum likelihood in order to simplify the problem. As in the case of the choice of λ svf , it might be also useful to estimate local LTI models on segments of training data to find out how the parameters are varying.

• Digital implementation of the CT filtering operations involved on the prefilter and auxiliary model: in order to obtain accurate filter outputs, the sampling frequency should be as high as possible.

• Numerical aspects: for recursive LS based methods, Potter's algorithm can be used to improve conditioning. For recursive IV based methods, this is not possible because the P matrix is non-symmetric. Therefore, the proposed solution is to use signal scaling. The advantage of this approach will be shown in Section 4.1.3.

• Estimator windup problem: the problem of having poor excitation in the data has been considered. For IV based methods, it is proposed to use conditional updating based on the condition number of the covariance matrix P (t k ), where P (t k ) is obtained from running in parallel RLSSVF with forgetting factor.

Problem formulation

According to the qualitative description of the circuit given in [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF], the system has two complex poles and one zero at the origin. Therefore, we consider the following CT LTV OE model

(p 2 + a 1 (t)p + a 2 (t))x(t) = b 0 (t)pu(t) (4.1a) y(t k ) = x(t k ) + e(t k ) (4.1b)
where the input u(t) and output y(t) are voltages. Given a data set, the aim is then to identify the parameters using the methods proposed previously in Chapters 1 and 2.

One of the assumption that we make to develop those approaches is that the parameters are slowly varying. Thus, from the many different scenarios available in [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF] 1 , we consider the data which satisfied that assumption. The name of the chosen data set is

MS_Harm_h3_N15640_RMS70_P2P700.mat (4.2)
which means the following:

• multisine excitation

• scheduling signal composed of one harmonic

• N = 15640 samples (in one period)

• RMS of 70 mV for the excitation

• scheduling signal with a peak-to-peak value of 700 mV

The data consists in 6 experiments, denoted here by Ei, with i = 1, . . . , 6. Each experiment has 3 periods, with N = 15640 samples in each period. In each experiment the input is different. Regarding the scheduling variable, there is one for the first 3 experiments, and another one for the last 3 experiments. We are using only the following experiments:

• E1 as estimation data

• E3 as validation data

• E4 as training data

Then, experiment E1 of the data set (4.2) will be denoted as (4.2)-E1.

The sampling frequency is 156250 Hz. Part of the input-output data is presented in Figure 4.2. The scheduling signal r(t) of the second period is shown in Figure 4.3. As we already mentioned, r(t) is not used in the identification. However it will be used to interpret the results.

To get a better insight of the system, we will also use the data set named MS_Const_N15640_RMS70_P2P700.mat

Recursive estimation of the electrical circuit

Our goal is to assess the performance of RLSSVF and RSRIVC. To get rough estimates of the hyperparameters λ svf and Q n , we first identify LTI models using the data set (4.3) corresponding to LTI systems. The estimation is performed with SRIVC with detrended data, i.e. the mean is removed. This procedure lead to λ svf = 10 7 rad/s. From [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF], we know that the parameter b 0 (t) is constant. Additionally we found out that also a 1 (t) is nearly constant. This means that in Q n the corresponding values for a 1 (t) and b 0 (t) are zero.

From the previous analysis, we know that there is only one value of Q n remaining to be estimated. This is done using the data set (4.2)-E1 corresponding to the LTV system. As part of the preprocessing, this data is also detrended. By trial and error we finally obtain Q n = diag([10 -15 0.5 10 -15 ]). Regarding t s,iv , from the convergence of RLSSVF, we choose t s,iv = 20 ms.

The RLSSVF estimates and RSRIVC estimates are compared in Figure 4.4. As expected, the results are very close due to the low noise level in the measurements. Notice that the variation of â2 (t) matches the shape of the scheduling signal r(t) displayed in Figure 4.3 for one period.

The frequency response functions for the frozen RSRIVC models corresponding to the second period are shown in Figure 4.5. Note that the shape is similar to the one reported in [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF]; as it is pointed out there, the evolution of the resonance frequency (peak of the frequency response function) matches the shape of the scheduling signal (see Figure 4.3).

To assess the quality of the estimates we use the benchmark criterion proposed in [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF], which corresponds to the RMS error computed on the second period of the simulated output signal, i.e.

RMSE y = 1 N -N tr 2N -1 k=N +Ntr y(t k ) -ŷ(t k , θ, u) (4.4)
where y is the measured output from the validation data. The simulated output ŷ is computed from the estimate θ and the measured input from the validation data u. The term N tr = 500 is used to eliminate the transient effects. Additionally, we compute the fit between the measured and simulated outputs, which is defined as follows

F = 100 × 1 - y(t k ) -ŷ(t k , θ, u) y(t k ) -mean(y(t k )) (4.5)
with t k varying as in (4.4). In [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF], for the analyzed data set, the proposed model structure, which is different than ours (see (4.1)), has 3 poles and one zero, with none of the parameters in the numerator being set to zero. Considering the data set (4.2)-E3 as validation data, we obtain after detrending the results that are shown in Table 4.1. From the values of RMSE y , we can see that RSRIVC is doing better than RLSSVF but still not as good as the result given in [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF]. Notice that in our case, apart from having a different model structure, we are only considering the data set (4.2)-E1 as estimation data. Nevertheless, from the fits we can conclude that the estimated models are a fair representation of the system dynamics.

For RSRIVC, a comparison between measured and simulated outputs is shown in Figure 4.6.

Illustration of the signal scaling approach

The advantage of using the signal scaling approach introduced in Section 3.5 is illustrated here. Signal scaling improves the conditioning of the proposed recursive algorithms, that is, it reduces the condition number of the P matrix.

The condition numbers of P corresponding to the previous estimations using RLSSVF and RSRIVC are shown in Figures 4.7(a) and (c), respectively. The large values obtained for κ 2 (P ) can be reduced through the signal scaling approach. We should keep in mind that by using signal scaling, the parameter estimates are not altered. The signal scaling approach requires the computation of the scaling matrix M s (see (3.53)), which is computed using the data set (4.2)-E4 after being detrended. The condition numbers of P with RLSSVF and RSRIVC with signal scaling decrease considerably, as we can see in Figure 4.7(b) and (d). Additionally, we provide κ 2 (P ) for the PotterSVF algorithm in Figure 4.7(e); as we can note, smaller condition numbers are obtained with the signal scaling approach.

Electronic throttle control system

In this section, the goal is to apply our proposed recursive estimation approaches to a simulated case study of an electronic throttle control (ETC). This corresponds to a closed-loop identification problem.

An ETC is a system used in automobiles to regulate air flow in the internal combustion engine. Basically, this is achieved by means of a DC (or step) motor which controls the position of the throttle valve plate. The ETC replaces the mechanical link between the vehicle acceleration pedal and the throttle valve in conventional automobiles. The benefits of the ETC, found both in gasoline and diesel engines, is to improve fuel economy, emissions and drivability. Additionally, the ETC is used for torque based engine control, where the desired torque is set by the pedal and the engine control system determines the necessary engine charge air and fuel to meet the torque requirement [Canudas de [START_REF] Canudas De Wit | Adaptive pulse control of electronic throttle[END_REF][START_REF] Pavković | Modeling and Control of Electronic Throttle Drive[END_REF][START_REF] Zhang | LPV Modeling and mixed constrained H 2 /H ∞ control of an electronic throttle[END_REF]].

An ETC system consists basically of a controller and an electronic throttle body (ETB), which is shown in Figure 4.8. Apart from the DC motor and throttle plate, there are a gearbox and two springs; to simplify the scheme, only one spring is shown in Figure 4.8. While the aim of the former is to reduce the speed coming from the motor, the aim of the latter is to return the plate into the so-called limp-home (LH) position in case of power supply failure [START_REF] Pavković | Modeling and Control of Electronic Throttle Drive[END_REF]. Non-linearities in the system come from the springs and also from friction effects. Additionally, there is the disturbance due to the air flow. On the other hand, the physical process parameters can vary, for instance, due to [START_REF] Pavković | Adaptive control of automotive electronic throttle[END_REF]]:

• Aging effect.

• Variations of external temperature which may produce changes in the armature resistance of the motor. The design of an ETC system is a challenging task, because of the ETB's nonlinearities and parameter variations. Among different approaches, a popular method for this application is sliding control (see e.g. [START_REF] Reichhartinger | Application of higher order sliding-mode concepts to a throttle actuator for gasoline engines[END_REF]). PID based controllers are also an option. In [START_REF] Deur | An electronic throttle control strategy including compensation of friction and limp-home effects[END_REF] for instance, a PID with compensation mechanisms for the nonlinearities has been developed. PIDs coupled with an adaptive mechanism have been also considered in [Canudas de [START_REF] Canudas De Wit | Adaptive pulse control of electronic throttle[END_REF][START_REF] Pavković | Adaptive control of automotive electronic throttle[END_REF].

In this chapter, we focus on the recursive estimation of the ETB parameters which could be used in an adaptive controller.

Electronic throttle body model

The present study is done based on the model provided by [START_REF] Bernardo | Synthesis and experimental validation of the novel LQ-NEMCSI adaptive strategy on an electronic throttle valve[END_REF], which represents an ETB for a gasoline engine. The ETB model is given by 4.9, where we can distinguish the electrical and mechanical submodels. The spring torque M s (t) is a piece-wise linear function. Here M s (t) is considered as defined in [di [START_REF] Bernardo | Synthesis and experimental validation of the novel LQ-NEMCSI adaptive strategy on an electronic throttle valve[END_REF]

di(t) dt = - R L i(t) - K v G r L ω(t) + 1 L u(t) (4.6a) ω(t) = dφ(t) dt (4.6b) J dω(t) dt = K t G r i(t) -M s (t) -M f (t) -M a (t
85 • π/180 rad J 1.56 • 10 -3 Kg•m 2 K b 4.98 • 10 -3 N•m•s/rad M c 8 • 10 -3 N•m R 1.4 Ω L 2.5 • 10 -3 H K v 0.02 V•s/rad K t 10.8856 • 10 -3 N•m/A G r 22.56 1
, i.e. M s (t) = h s (φ(t)) =              h s3 (φ(t)) if φ ∈ [φ min ; φ LH -∆φ/2[ 0 if φ ∈ [φ LH -∆φ/2; φ LH + ∆φ/2] h s1 (φ(t)) if φ ∈ ]φ LH + ∆φ/2; φ 12 [ h s2 (φ(t)) if φ ∈ [φ 12 ; φ max ] (4.7) with h s1 (φ(t)) = K s1 φ(t) -φ LH + ∆φ 2 + M open (4.8a) h s2 (φ(t)) = K s2 (φ(t) -φ 12 ) + h s1 (φ 12 ) (4.8b) h s3 (φ(t)) = -K s3 φ LH - ∆φ 2 -φ(t) -M close (4.8c)
The values of the coefficients in (4.7) and (4.8) are given in Table 4.3. A plot of M s (t) is presented in Figure 4.10.

The literature about friction modeling is extensive (see e.g. [START_REF] Olsson | Friction models and friction compensation[END_REF]] and references therein). To represent friction in an ETB, both static and dynamical models have been used. In [di [START_REF] Bernardo | Synthesis and experimental validation of the novel LQ-NEMCSI adaptive strategy on an electronic throttle valve[END_REF], the static model consists of Coulomb, Stribeck and viscous friction. Regarding a dynamic model, the so-called LuGre model is considered in [Canudas de [START_REF] Canudas De Wit | Adaptive pulse control of electronic throttle[END_REF]. In our case, the model proposed in [di [START_REF] Bernardo | Synthesis and experimental validation of the novel LQ-NEMCSI adaptive strategy on an electronic throttle valve[END_REF] is simplified, such that the Stribeck effect is neglected. Then, M f (t) is given by

M f (t) = h f (ω(t))
= K b ω(t) + M c sign(ω(t)) (4.9)

with the coefficients defined in Table 4.3.

The aerodynamic torque M a (t) depends e.g. on the throttle plate angle φ(t) among others [Canudas de [START_REF] Canudas De Wit | Adaptive pulse control of electronic throttle[END_REF]]. Here we consider that it is zero.

Problem formulation

We define now the true system for the simulations based on the model (4.6). As it was mentioned before, the system is subject to parameter variations. Some changes can be very slow due to aging effects. However, this is not the case for the armature resistance, which can differ considerably within a single cycle of engine operation, due to temperature variations in the DC windings [Pavković and Deur, 2011, p. 97]. According to [Pavković and Deur, 2011, p. 98], the armature resistance does not exceed ±50% of the nominal value. The same variation is assumed here, i.e.

0.5 • R ≤ R o (t) ≤ 1.5 • R (4.10)
with the nominal value R specified in Table 4.3. The rest of the true parameters, denoted also with a superscript {} o , are defined according to Table 4.3. Since the electrical dynamics is much faster than the mechanical dynamics, the dynamics of the current is usually neglected. Such an assumption is also considered here both for the true system and the model to be presented next. In this way we can better evaluate the tracking of the parameters. Then, the true system S is given by The system operates in closed loop, and the closed-loop configuration is shown in Figure 4.11. The true plant G o is defined by the first three equations in (4.11). F φ is a filter to smooth the set-point and C is a PID controller. The former is described by F φ (q -1 ) = 0.009901 + 0.009901q -1 1 -0.9802q -1 (4.12) and the latter by

i(t) = 1 R o (t) u(t) - K o v G o r R o (
C : C c (q -1 ) = k p + k i 1 -q -1 + k d (1 -q -1 ) (4.13)
with k p = 2.85, k i = 6.06 • 10 -3 , k d = -36.96.

Due to the complexity of the system, the model identification of a throttle valve is done considering more than one experiment (see e.g. [START_REF] Isermann | Fault diagnosis of machines via parameter estimation and knowledge processing -Tutorial paper[END_REF]Münchhof, 2010, Pavković and[START_REF] Pavković | Modeling and Control of Electronic Throttle Drive[END_REF]). In this study, assuming that the parameters of the mechanical submodel have been previously identified, we focus on the recursive estimation of the electrical submodel which is given by with e(t k ) a zero-mean DT Gaussian noise sequence. Note that both the electrical and mechanical subsystem have the parameter G o r , and it is assumed to be known. The parameter vector to be estimated is i.e. we need to identify the linear function h s1 (φ(t)) defined in (4.8a). Therefore, the non-linear effects are only due to friction. Notice that the closed-loop system is excited only through the set point φ d (t), i.e. no external signal2 is considered, since then it is easier to satisfy (4.16) by a proper choice of φ d (t). Then, the identification problem is to recursively estimate the time-varying parameter vector (4.15) that characterizes the model (4.14), based on the sequences

i(t) = 1 R(t) u(t) - K v G o r R(
θ(t) = θ 1 (t) θ 2 (t) T = 1 R(t)
{φ d (t k ), u(t k ), i(t k ), φ m (t k )} N ′ k=1
where N ′ is the number of samples which increases by one with every recursion.

Recursive estimation of the electrical submodel

Let us formulate the solution with CLRSRIVC. Firstly, note that (4.14) is not a CT LTV OE model, which is the type of model considered in the development of the identification Figure 4.12: Closed-loop auxiliary model. methods in Chapter 2. However, the proposed methods can be adapted to handle this situation.

From a heuristic approach we found out that a suitable prefilter for the identification of the electrical model is the prefilter corresponding to the mechanical system, which can be written as follows Note that if we discard the last term in (4.17), which is the source of the non-linearity, we obtain a multi-input single-output linear time-varying (MISO LTV) model with 2 inputs. The prefilter corresponding to this MISO LTV model is

p 2 + K o b J o p + K o s1 J o φ(t) = K o t G o r J o i(t) - m o s1 J o - M o c J o sign
F (p) = K o s1 /J o p 2 + (K o b /J o )p + K o s1 /J o (4.19)
Since we assume that the parameters of the mechanical subsystem are known and constant, then the prefilter (4.19) is also known and constant. Now we rewrite (4.14) as a filtered linear regression. By assuming slow parameter variations, we can neglect the non-commutativity issue and then apply the prefilter (4.19) to (4.14a) obtaining i f (t) = θ 1 (t)u f (t) -θ 2 (t)φ The closed-loop auxiliary model to compute the instruments is shown in Figure 4.12, where the estimated plant Ĝ is given by î Chapter 5

(t) = 1 R(t) û(t) - Kv G o r R (
Conclusions and future work

Conclusions

In this thesis we have addressed the on-line identification of CT LTV models. In particular we focused on CT output-error type models whose parameters vary slowly in time.

In Chapter 1, we have first shown how linear filter methods are able to estimate CT LTI models. A simple approach that has been recalled is the LSSVF method which gives a closed form solution. That is what makes this least squares based method attractive. However, if the output measurement noise is not negligible, the LSSVF estimates will be biased. For the identification of CT LTV models, the recursive counterpart of LSSVF, denoted RLSSVF, can be used. In LTV model identification, analogously to the LTI case, the presence of measurement noise deteriorates the RLSSVF estimates, and the degradation of the results might not be acceptable for certain applications.

In Chapter 2, the contribution has been to develop recursive instrumental variable based methods, both for open-loop and closed-loop identification of CT LTV models. As RLSSVF, these are also direct identification methods based on the linear filter approach, but provide better estimates when the output measurement noise is not negligible. To track the time-varying parameter variations, we proposed the forgetting factor and the random walk model approaches. Nevertheless, we have focused on the latter since it is more flexible, meaning that it is suitable for cases in which the parameters vary at different rates. Furthermore, it should be more robust in cases of poor excitation. Numerical examples including Monte Carlo simulations have been considered to show that IV methods are advantageous over RLSSVF in terms of delivering parameter estimates with smaller tracking errors. One of the numerical examples and some of the algorithms developed in this chapter, have been included in the CONTSID toolbox.

In Chapter 3 we have proposed solutions to the practical aspects and implementation issues related to the use of the recursive IV methods that have been developed. The main contributions are:

• Development of a method to estimate the normalized covariance matrix of the Kalman filter Q n . This is one of the hyperparameters of the developed identification methods that has to be chosen by the user. The tracking ability and noise sensitivity of Kalman filter based algorithms depend on Q n . If the parameters of the system are not varying at the same speed, different values for the diagonal elements of Q n have to be chosen. Finding suitable values for the covariance matrix might be time-consuming and not straightforward. We have proposed an off-line method to estimate Q n from a given training data set. The approach couples the maximum likelihood method with one of the developed IV based methods. The effectiveness of this approach has been shown in a numerical example.

• Proposition of guidelines to choose the hyperparameters of the developed identification methods. Guidelines were given for the following hyperparameters, which along with Q n , are the most important:

-Cut-off frequency of the SVF (λ svf ). This is an hyperparameter of RLSSVF which is used to initialize the IV based methods. The impact of λ svf in the estimates have been discussed. In particular, we showed through a numerical example how overshoots of the estimates in the transient mode are related to λ svf . Approaches to reduce these overshoots have been also presented.

-Time instant when the IV based method starts to operate (t s,iv ). The value of t s,iv should be the time instant at which the RLSSVF has converged. Several measures that can be used to choose t s,iv were discussed.

• Proposition of an anti-windup approach for IV based methods. When a system to be identified is poorly excited, anti-windup methods are needed to prevent the parameter estimates to blow up. For the IV based methods, the solution that we have proposed is to use conditional updating based on the condition number of the parameter covariance matrix P , where P is obtained from the RLSSVF with forgetting factor algorithm running in parallel. The motivation to use the RLSSVF with forgetting factor method is that it should be in general more sensitive than other algorithms to poor excitation.

• Development of a signal scaling approach to improve conditioning of IV based methods. Numerical aspects can play an important role when recursive algorithms are run in a microcontroller. To improve the conditioning of the identification problem and thus to avoid numerical issues, we have proposed to use a signal scaling approach. The effectiveness of this method has been demonstrated with a benchmark considering real data from an electronic band pass filter.

Finally, in Chapter 4, a simulation case study corresponding to the estimation of an electronic throttle valve model is considered. The identification problem is challenging because of non-linear effects due to friction. Moreover the system operates in closedloop. We have illustrated how to adapt one of the developed IV methods to the recursive estimation of the physical time-varying parameters.

Future work

In this thesis we have developed recursive instrumental variable methods for CT models under the assumption that the parameters are slowly varying in time. This assumption allows us to neglect the non-commutativity problem when linear filter method are applied. Additionally this assumption is used to guarantee stability of the auxiliary model. The goal of a future work can be to adapt the developed approaches to scenarios with fast or abrupt parameter variations. The problem is not trivial because of these non-commutativity and stability issues.

For fast parameter variations, a solution might be to use a parallel estimation scheme, were several recursive algorithms running in parallel compete with each other. These algorithms could differ in the assumed representation of the time-varying parameters which can be described by different deterministic and/or stochastic models. For the case of abrupt parameter variations, parallel estimation is also an option.

One of our algorithms, developed under the assumption of a having a CT LTV system, has been adapted to a CT non-linear time-varying system. Another possible direction of research would be to develop algorithms specially dedicated to this type of systems. Notice that (C.9) involves many multiplications with the diagonal matrix M s , which means that several multiplications with zeros are performed. A more efficient scalar implementation could be developed.

Résumé

Les travaux présentés dans ce mémoire traitent de l'identification des systèmes dynamiques représentés sous la forme de modèles linéaires continus à paramètres variant lentement au cours du temps. La complexité du problème d'identification provient d'une part du caractère inconnu de la loi de variation des paramètres et d'autre part de la présence de bruits de nature inconnue sur les signaux mesurés.

Les solutions proposées s'appuient sur une combinaison judicieuse du filtre de Kalman en supposant que les variations des paramètres peuvent être représentées sous la forme d'une marche aléatoire et de la méthode de la variable instrumentale qui présente l'avantage d'être robuste vis à vis de la nature des bruits de mesure.

Les algorithmes de type récursif sont développés dans un contexte d'identification en boucle ouverte et en boucle fermée. Les différentes variantes se distinguent par la manière dont est construit la variable instrumentale. Inspirée de la solution développée pour les systèmes linéaires à temps invariant, une construction adaptative de la variable instrumentale est suggérée pour pouvoir suivre au mieux l'évolution des paramètres. Les performances des méthodes développées sont évaluées à l'aide de simulations de Monte Carlo et montrent la suprématie des solutions proposées s'appuyant sur la variable instrumentale par rapport celles plus classiques des moindres carrés récursifs.

Les aspects pratiques et d'implantation numérique sont d'une importance capitale pour obtenir de bonnes performances lorsque ces estimateurs sont embarqués. Ces aspects sont étudiés en détails et plusieurs solutions sont proposées non seulement pour robustifier les estimateurs vis à vis du choix des hyper-paramètres mais également vis à vis de leur implantation numérique. Les algorithmes développés sont venus enrichir les fonctions de la boîte à outils CONTSID pour Matlab.

Enfin, les estimateurs développés sont exploités pour faire le suivi de paramètres de deux systèmes physiques : un benchmark disponible dans la littérature constitué d'un filtre électronique passe-bande et une vanne papillon équipant les moteurs de voiture. Les deux applications montrent le potentiel des approches proposées pour faire le suivi de paramètres physiques variant lentement dans le temps.

Mots-clés : modèle à temps continu, modèle linéaire à paramètres variables, estimateur de la variable instrumentale, filtre de Kalman

Figure 1 :

 1 Figure 1: Entités de base pour l'identification d'un système: jeu de données, structure de modèles et méthode d'identification.

Figure 2 :

 2 Figure 2: Approches indirectes et directes pour l'identification d'un modèle CT.

1. 3 .

 3 Figure 1.3: Indirect and direct recursive identification of a physical parameter vector θp (t k ).

Figure 1

 1 Figure 1.4: Part of the input-output data.

Figure 1 Figure 1 Figure 1

 111 Figure1.6: Bode plots for the frozen systems obtained every 250 s. The bode plots for t 1 = 0 and t N ′ = 9000 s, matching in this case, are highlighted.

Figure 1

 1 Figure 1.9: True CT parameters and CT estimates from the indirect RLS and direct RLSSVF methods.

A

  . Padilla, H. Garnier, and M. Gilson. Version 7.0 of the CONTSID toolbox, 17th IFAC Symposium on System Identification, Beijing China, October 2015. A. Padilla, H. Garnier, P. C. Young, J. I. Yuz, Real-time identification of continuoustime slowly linear time-varying systems, 55th IEEE Conference on Decision and Control, Las Vegas US, December 2016. A. Padilla, H. Garnier, P. C. Young, J. I. Yuz, Recursive online IV method for identification of continuous-time slowly time-varying models in closed loop, IFAC World Congress, Toulouse France, July 2017.Chapter 2

Figure 2

 2 Figure 2.2: Diagram of Kalman filter based RSRIVC.

  ) and ζ f (t k ) are computed as in the recursive LTI case (see Section 2.2.3). The methods are called Kalman filter based RIVSVF and Kalman filter based RSRIVC algorithms. A diagram of the Kalman filter based RSRIVC method is presented in Figure 2.2; note that the auxiliary model and adaptive filter are shown in their operator form.

Figure 2 . 3 :

 23 Figure 2.3: True parameters, RLSSVF estimates and RSRIVC estimates for experiment 1.

Figure 2

 2 Figure 2.4: Average of the parameter relative errors for RLSSVF and RSRIVC for a Monte Carlo simulation with 100 runs.

Figure 2 Figure 2 . 6 :Figure 2

 2262 Figure 2.5: Absolute value of lag errors for RLSSVF and RSRIVC w.r.t. experiment 1.

  where a o 1 (t) varies slowly between 5 and 45 in a linear fashion, b o 0 (t) remains constant at 200 and a o 2 (t) varies slowly as follows, a o 2 (t) = 160 -90 cos(2πt/1000)

Figure 2 Figure 2 . 9 :Figure 2 Figure 2

 22922 Figure 2.8: Magnitude Bode diagram for some frozen systems obtained every 100 s. The system responses for t = 0 and t = 1000 s are highlighted.

Figure 2 Figure 2 Figure 2

 222 Figure 2.12: True parameters, RIVSVF estimates (a) and RSRIVC estimates (b).

Figure 2 .

 2 Figure 2.15: Closed-loop system.

Figure 2 .

 2 Figure 2.16: Auxiliary model.

Figure 2 Figure 2

 22 Figure2.17: True parameters and RLSSVF estimates for the noise-free case and noisy case. (The estimates for the deterministic case are matching the true values).

Figure 2

 2 Figure 2.19: Average of the relative errors for RLSSVF for a Monte Carlo simulation with 100 runs.

Figure 3

 3 Figure 3.1: Diagram of RSRIVC with issues (Ii) regarding practical usage and implementation (see text for the description of the issues).

Figure 3 . 4 :Figure 3

 343 Figure 3.3: Part of the input-output data.

Figure 3

 3 Figure 3.6: Part of the input-output data.

Figure 3 Figure 3 . 8 :

 338 Figure 3.7: True parameters and RLSSVF estimates with and without initial condition handling.

Figure 3 . 9 :

 39 Figure 3.9: True parameters and RLSSVF estimates using t s,ls = 2 s.

Figure 3

 3 Figure 3.10: Part of the input-output data.

  includes experiments with constant scheduling variables with different values.

Figure 4 . 2 :Figure 4 . 3 :Figure 4

 42434 Figure 4.2: Part of the input-output data for the identification of the circuit model.

Figure 4

 4 Figure 4.5: Frequency response functions of the frozen models obtain with RSRIVC.

Figure 4 . 6 :Figure 4 . 7 :

 4647 Figure 4.6: Comparison between measured and simulated outputs for RSRIVC (the measured and simulated outputs are nearly matching).

Figure 4 . 8 :

 48 Figure 4.8: Electronic throttle body (ETB).

Figure 4 . 9 :Figure 4

 494 Figure 4.9: Block diagram of the ETB model.

Figure 4 .

 4 Figure 4.11: Closed-loop block diagram of the ETB.

  avoid non-linearities coming from the spring, the estimation data satisfies the condition φ ∈ ]φ LH + ∆φ/2; φ 12 [ = ]13.14; 85[ (deg) (4.16)

Figure 4 Figure 4 Figure 4

 444 Figure 4.14: True parameters, RLSSVF and CLRSRIVC estimates for a single experiment.

Table 2 .

 2 1: Mean squared value of the parameter tracking error.

	Method	MSE[ θ(t k )]
	RLSSVF	1.152
	RIVSVF	0.814
	RSRIVC	0.800

Table 2 .

 2 2: Mean squared of the lag error and mean squared of the noise error w.r.t.

	experiment 1.
	Method	MSE l [ θ(t k )] MSE 1 n [ θ(t k )]
	RLSSVF	0.8036	0.2578
	RIVSVF	0.4114	0.2670
	RSRIVC	0.3998	0.2672

Table 3 .

 3 1: Overshoots for the estimates â1 and b0 .

	Method & settings	overshoots (%) for:
		â1	b0
	RLSSVF without i.c.	6561
	RLSSVF with i.c.	110	0.005
	RLSSVF with t s,ls = 2 s 0.0052	0.0193

Table 4 .

 4 1: Performance indexes using validation data.

	Method	RMSE y (mV)	F (%)
	[Lataire et al., 2015]	0.1156	not reported
	RLSSVF	0.6909	95.7
	RSRIVC	0.3781	97.6

Table 4 .

 4 2: Description of ETB model parameters.

	Symbol Description
	φ	Angular position of the throttle plate
	ω	Angular velocity of the throttle plate
	i	Current
	u	Voltage source of the armature
	L	Inductance of the armature coil
	R	Resistance of the armature coil
	K v	Velocity coefficient
	K t	Torque coefficient
	J	Moment of inertia
	G r	Transmission ratio due to the gear
	M s	Spring torque
	M f	Friction torque
	M a	Aerodynamic torque due to the air flow

Table 4 .

 4 3: ETB physical parameter values coming from [di[START_REF] Bernardo | Synthesis and experimental validation of the novel LQ-NEMCSI adaptive strategy on an electronic throttle valve[END_REF].

	Symbol	Value	Unit
	K s1	1.83 • 10 -3 • 180/π	N•m/rad
	K s2	22.90 • 10 -3 • 180/π N•m/rad
	K s3	47.15 • 10 -3 • 180/π N•m/rad
	M close	28.4162 • 10 -3	N•m
	M open	15.8202 • 10 -3	N•m
	φ min	8 • π/180	rad
	φ max	90 • π/180	rad
	φ LH	12.58 • π/180	rad
	∆φ	1.12 • π/180	rad
	φ 12

  (t k ) = φ(t k ) + e o (t k ) (4.11d)with R o (t) varying linearly in the range (4.10). Regarding the signal e o (t k ), it is a zero-mean DT Gaussian noise sequence with variance 10 -3 (rad/s) 2 . The torques M s (t) and M f (t) are also defined considering the values provided in Table4.3 as the true parameters.

			t)	ω(t)	(4.11a)
	ω(t) =	dφ(t) dt	(4.11b)
	J o dω(t) dt	= K o t G o r i(t) -M s (t) -M f (t)	(4.11c)
	φ m

  Ms (t) and Mf (t) defined with the true parameters. The other chosen hyperparameters are: t s,iv = 5 s; Q n = diag([10 -6 5 • 10 -5 ]).For comparative purposes we also consider RLSSVF with the SVF

			t)	ω(t)	(4.21a)
	ω(t) = J o d ω(t) dt = K o d φ(t) dt t G o r î(t) -Ms (t) -Mf (t)	(4.21b) (4.21c)
	with the torques F (p) =	λ svf p + λ svf	2		(4.22)

  Then, the scaled version of the Kalman filter based Potter's algorithm is given by Prediction step:θs (t k |t k-1 ) = θs (t k-1 ) (C.9a) (t k ) -ϕ f T (t k )M s θs (t k |t k-1 ) (C.9d) g(t k ) = S T s (t k |t k-1 )M s ϕ f (t k ) (C.9e) β(t k ) = 1 + g T (t k )g(t k ) (C.9f) α(t k ) = β(t k ) + β(t k ) (C.9g) L(t k ) = M s S(t k |t k-1 )S T (t k |t k-1 )M s ϕ f (t k )

	S T s (t k |t k-1 ) 0	= T	S T s (t k-1 ) Q n,s T /2	(C.9b)
	Correction step:
	θs (t k ) = θs (t k |t k-1 ) + M -1 s L(t k )ε(t k )	(C.9c)
	ε(t k ) = y	(na) f
				β(t k )	(C.9h)
	S

s (t k ) = S s (t k |t k-1 ) I d -g(t k )g T (t k ) α(t k ) (C.9i)

The static gain can be made unity.

A CT LTV OE model set will be considered in Chapter 2.

www.cran.univ-lorraine.fr/contsid/

An analysis about signal scaling and its impact in the identification problem is presented in Section 3.5.

Unless otherwise stated, we consider M θ (q -1 ) = q -1 both in the previous and the following examples.

The condition number of a matrix A, denoted by κ2(A), is defined in Section 3.5.1.

The data can be downloaded from www.kth.se/social/group/system-identificatio/page/17th-ifac-symposium-on-system-identifica/

The external signal corresponds to r1(t) in Figure2.15.

Acknowledgements

Chapter 4

Applications

Here we present two real-life problems to which we apply the approaches developed in Chapters 2 and 3. The first one is an electronic bandpass filter proposed as a benchmark both for LPV and LTV model identification [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF]. The second real life based problem is an electronic throttle control (ETC) system that we find in automobiles nowadays. In the last decades, several studies have been carried out in order to improve the performance of this device (see e.g. [START_REF] Pavković | Modeling and Control of Electronic Throttle Drive[END_REF]).

Benchmark data from an electronic bandpass filter

In [START_REF] Lataire | Benchmark data on a linear time-and parameter-varying system[END_REF], real data collected from a second order bandpass filter has been proposed as a benchmark both for LPV model and LTV model identification. It deals with an open-loop estimation problem. In Figure 4.1, a diagram of the electronic circuit is shown. A parallel connection of an n-type J-FET transistor and a 470 kΩ resistor constitute a variable resistance, which produces the parameter variations. The scheduling variable r(t) is the gate-source voltage of the transistor, and the input u(t) and output y(t) are the voltages indicated in Figure 4.1. The noise level in the measurements is very low, with an SNR of more than 60 dB. Note that some small non-linear effects are also present.

Several scenarios, consisting on different excitations and scheduling signals are provided. To assess our approaches, we consider one scenario with smooth parameter variations. Since LTV models are identified, the measured scheduling signal is not used. where λ svf is defined in terms of the bandwidth ω b of (4.17), which is 12.4 (rad/s); the chosen value is λ svf = 5.0 • ω b . For Q n we take the same matrix as for CLRSRIVC. We present next some results for a single experiment run out of the 100 that are considered later in a Monte Carlo simulation analysis. The measurements that are used for the model estimation are shown in Figure 4.13. The results for RLSSVF and CLRSRIVC are presented in Figures 4.14 (a) and (b), respectively. We can clearly see that while for RLSSVF there is a significant tracking error, for CLRSRIVC the estimates are very close to the true values. From θ(t) it is straightforward to compute the physical estimates R(t) and Kv (see (4.15)). A comparison between the true physical parameters and the CLRSRIVC estimates is shown in Figure 4.15.

The Monte Carlo simulation is performed with 100 experiments. The results, presented in Figure 4.16, confirm the previous analysis for a single experiment.

Appendix A

Digital implementation of CT filtering operations

The digital implementation issues of the CT filtering operations are well-known in CT model identification. They should be treated in an appropriate way since errors generated by the digital implementation can have an impact on the quality of the estimated model.

The CT filtering operation is represented by a LTV state-space system

and J(t) of appropiate dimensions. Depending on the intersample behaviour of the input u(t), two different discretizations are considered in this work.

A.1 Zero-order-hold assumption

The zero-order-hold assumption is considered if the intersample behaviour of the input u(t) is piecewise constant. In the LTV scenario, the zero-order-hold assumption applies also to the system time-varying matrices F (t), G(t), H(t) and J(t), i.e. the following is satisfied for the input

and for the time-varying matrices we have

That is, between samples the matrices are constant. As pointed out in [START_REF] Tóth | Discretisation of linear parametervarying state-space representations[END_REF] in the context of LPV systems, this assumption of having constant matrices is usually not realistic. The discrete-time version of (A.1) is then defined as in the LTI case (see [Åstrom and Wittenmark, 1996, p. 34]) and it is given by

where

Ts 0 e F ((k-1)Ts)s ds G (A.5b)

These matrices F d (t k ) and G d (t k ) can be computed from

with 'exp' corresponding to the matrix exponential that can be computed in MATLAB with expm.

A.2 Trapezoidal approximation

If the intersample behaviour for the input u(t) is different than piecewise constant, the trapezoidal approximation may be used assuming that the system time-varying matrices F (t), G(t), H(t) and J(t) remain constant between samples, i.e. (A.3) is satisfied.

The discrete-time version of (A.1) is then defined as in the LTI case (see [Franklin et al., 1998, p. 200]) and it given by (A.4), where

Stability of LTV systems

The following material is extracted mainly from [Khalil, 2002, pp. 156-158]. Let us consider the LTV system ẋ(t) = F (t)x(t) (B.1) with state x(t) ∈ R nx and system matrix F (t). The solution of (B.1) is given by

where Φ(t, t 0 ) is the state transition matrix. The following theorem defines stability in terms of the state transition matrix.

Theorem 1 The equilibrium point x = 0 of (B.1) is uniformly exponentially stable if and only if the state transition matrix satisfies the inequality

for some positive constants k and λ.

To test stability in practice, instead of using Theorem 1, it is more convenient to work with the following theorem.

Theorem 2 Let,

• x = 0 be the exponentially stable equilibrium point of (B.1).

• Q(t) be a continuous, bounded, positive definite, symmetric matrix.

Then, there is a continuously differentiable, bounded, positive definite, symmetric matrix

x is a Lyapunov function for the system that satisfies the conditions of Theorem 1.

Equation (B.

3) corresponds to a linear matrix inequality that can be solved numerically.

Appendix C

Scaled version of Potter's algorithm

For simplicity of exposition, the time index is omitted.

), for the former we use the superindex {•} + . Then, the a-posteriori covariance matrix is written as follows

Replacing (3.65b) in (C.1) we obtain

Using the square root factorization P s = S s S T s yields

then (C.3) can be written in the form of (3.69), i.e.

Since (C.5) has the same form than (C.3), we know that

, and U can be written as follows

where α is given by (3.72) but with β defined in (C.6). Finally, considering (C.4a), the square root of P s can be recursively computed as follows

Abstract

The work presented in this thesis deals with the identification of dynamic systems represented through continuous-time linear models with slowly time-varying parameters. The complexity of the identification problem comes on the one hand from the unknown character of the parameter variations and on the other hand from the presence of noises of unknown nature on the measured signals.

The proposed solutions rely on a judicious combination of the Kalman filter assuming that the variations of the parameters can be represented in the form of a random walk, and the method of the instrumental variable which has the advantage of being robust with respect to the nature of the measurement noises.

The recursive algorithms are developed in an open-loop and closed-loop identification setting. The different variants are distinguished by the way in which the instrumental variable is built. Inspired by the solution developed for time-invariant linear systems, an adaptive construction of the instrumental variable is suggested in order to be able to follow the evolution of the parameters as well as possible. The performance of the developed methods are evaluated using Monte Carlo simulations and show the supremacy of the proposed solutions based on the instrumental variable compared with the more classical least squares based approaches.

The practical aspects and implementation issues are of paramount importance to obtain a good performance when these estimators are used. These aspects are studied in detail and several solutions are proposed not only to robustify the estimators with respect to the choice of hyperparameters but also with respect to their numerical implementation. The algorithms developed have enhanced the functions of the CONTSID toolbox for Matlab.

Finally, the developed estimators are considered in order to track parameters of two physical systems: a benchmark available in the literature consisting of a bandpass electronic filter and a throttle valve equipping the car engines. Both applications show the potential of the proposed approaches to track physical parameters that vary slowly over time.