Michel Juillard
email: michel.juillard@ens.fr

Sébastien Villemot
email: sebastien.villemot@ens.fr

Multi-country real business cycle models: Accuracy tests and test bench

This paper describes the methodology used to compare the results of different solution algorithms for a multi-country real business cycle model. It covers in detail the structure of the model, the choice of values for the parameters, the accuracy tests used in the comparison, and the computer program specifically developed for performing the tests.

Introduction

In order to compare complex algorithms such as the ones necessary to compute the solution of nonlinear dynamic stochastic general equilibrium models, one needs to define clearly the dimensions along which the results will be compared. Furthermore, it is also important to compare the results of the different algorithms using the same implementation of the tests so that one can exclude differences in the administration of the testing procedure. For this reason, separately from the implementation of each algorithm by the participants to the comparison project, we have developed a computer program that takes as input the solution provided by each group and subjects it to an identical set of tests.

Following previous literature (i.e. [START_REF] Aruoba | Comparing solution methods of dynamic equilibrium economies[END_REF][START_REF] Den Haan | Accuracy in Simulations[END_REF][START_REF] Heer | Dynamic general equilibrium modelling, computational methods and applications[END_REF][START_REF] Judd | Projection methods for solving aggregate growth models[END_REF], we focus on evaluating errors of approximation: the errors that one obtains by plugging back the 1 approximated solution in the original problem. We implement this approach in two different setups: the first test computes errors of approximation at different points in the state space located on a sphere around the deterministic steady state of the problem. This test informs us on the ability of a given solution method to provide accurate approximations of the solution further and further away from a point "easy" to solve for: the deterministic steady state. Of course, ordering points by their distance to the steady state does not indicate how likely these points are in simulated data. This remark motivates the second test: measuring the errors of approximation along a given simulation path. More precisely, each solution is simulated starting from the same initial point, the deterministic steady state, and subjected to the same set of random shocks. In the literature, the test proposed by Den [START_REF] Den Haan | Accuracy in Simulations[END_REF] is often used to test whether prediction errors implied by the numerical solutions are orthogonal to the elements of information set. Finally, we didn't retain it in the testing procedure used for this exercise, for reasons that we present below. All approaches provide an approximated solution to the first-order conditions of the social planner optimization problem. In the first section, we provide a formal presentation of the model and derive the first-order conditions of the social planner problem. The detailed specifications and the parameterization of the 30 models are given in section 2. The approximation errors are defined in section 3 along with the specification of the tests retained in the comparison. In section 4, we describe the computer program that serves as the test bench. It is this software that is used in the comparison paper by [START_REF] Kollmann | Comparison of numerical solutions to a suite of multi-country real business cycle models[END_REF].

1 Description of the model economies

General setup

We consider a multi-country RBC model with complete asset markets. We denote by N the number of countries. A single homogenous good is produced, traded and consumed across countries.

Production of country j ∈ {1, . . . , N} at date t is equal to a j t f j (k j t , j t ), where f j is the production function, k j t is the beginning-of-period capital stock, j t is hours worked, and a j t is productivity level. The law of motion of capital is:

k j t+1 = (1 -δ)k j t + i j t (1)
where i j t is investment and δ is the depreciation rate of capital. The law of motion of productivity is:

ln a j t = ρ ln a j t-1 + σ(e j t + e t ) ( 2 )
where e j t is a country specific shock and e t is a worldwide shock. Shocks are assumed to be independent and identically distributed Gaussian variables with zero mean and unit variance. Parameters ρ and σ determine the autocorrelation and variance of the logarithm of the productivity level, respectively.

There is an adjustment cost on the capital stock:

Γ j t = φ 2 k j t i j t k j t -δ 2 (3)
where parameter φ governs the intensity of the friction. The value of this parameter has an impact on the dimension of the state space. When φ > 0,

all the country-level capital stocks are state variables. In contrast, when φ = 0 (i.e. when there is no adjustment cost), only the world capital stock matters, since in that case capital is perfectly mobile. In this project, we are only interested in specifications where the stock of capital of each country matters (φ > 0).

Each country has a representative agent, whose instantaneous utility is

u j (c j t , j t )
, where u j is the utility function and c j t is consumption.

The aggregate world resource constraint is:

N j=1 c j t + i j t -δ k j t ≤ N j=1 a j t f j (k j t , j t ) -Γ j t (4)
It is a well known result that the decentralized market equilibrium with complete asset markets is equivalent to the solution of a corresponding social planner problem, where each country has a weight τ j (Negishi weight) in the planner's objective, the weights depending on initial endowments.

We therefore formulate the problem as a social planner problem. Since we want to evaluate algorithms used for general equilibrium problems, solution algorithms are required to use Euler equations together with dynamic and static constraints to compute the decision rules; in particular, this rules out algorithms that compute the value function of the social planner.

The problem of the social planner is: max

{c j t ,i j t ,k j t+1 , j t } j=1,...,N t=0,...,+∞ E 0 N j=1 τ j +∞ t=0 β t u j (c j t , j t ) (5)
subject to (1), ( 3) and ( 4), where β is the subjective time discount factor.

We call λ t the Lagrange multiplier of the aggregate resource constraint, which is binding in equilibrium. The first-order conditions are:

τ j u j c (c j t , j t ) = λ t (6) τ j u j (c j t , j t ) = -λ t a j t f j (k j t , j t ) ( 7 ) λ t 1 + φ i j t k j t -δ = β E t λ t+1 1 + a j t+1 f j k (k j t+1 , j t+1 ) ( 8 ) +φ 1 -δ + i j t+1 k j t+1 - 1 2 i j t+1 k j t+1 -δ i j t+1 k j t+1 -δ k j t+1 = (1 -δ)k j t + i j t (9) N j=1 c j t + i j t -δk j t = N j=1 ⎡ ⎣ a j t f j (k j t , j t ) - φ 2 k j t i j t k j t -δ 2 ⎤ ⎦ (10) ln a j t = ρ ln a j t-1 + σ e t + e j t ( 11
)
where u j c (u j ) is the derivative of u j with respect to c ( ), and f j k (f j ) is the derivative of f j with respect to k ( ).

Specifications of utility and production functions

Four specifications are used for the utility function:

1. A constant relative risk aversion (CRRA) function, incorporating only consumption (no labor):

u j (c j t , j t ) = (c j t ) 1-1 γ j 1 -1 γ j (12)
where γ j is the intertemporal elasticity of substitution.

2. A utility function separable in consumption and labor:

u j (c j t , j t ) = (c j t ) 1-1 γ j 1 -1 γ j -b j ( j t ) 1+ 1 η j 1 + 1 η j (13)
where η j is the Frisch elasticity of labor supply, and b j governs the relative weight of consumption and leisure.

3. A Cobb-Douglas utility function (hence with unit elasticity of substitution between consumption and leisure):

u j (c j t , j t ) = γ j γ j -1 (c j t ) ψ j L e -j t (1-ψ j ) 1-1 γ j (14)
where L e is the labor endowment of the representative agent in each country, and ψ j governs the relative weight of consumption and leisure.

A constant elasticity of substitution (CES) function embedded inside

a CRRA function:

u j (c j t , j t ) = γ j γ j -1 (c j t ) 1-1 χ j + b j L e -j t 1-1 χ j 1-1 γ j 1-1 χ j (15)
where χ j is the elasticity of substitution between consumption and leisure.

Three specifications are used for the production function:

1. A production function with decreasing returns to scale, using only capital:

f j (k j t , j t ) = A(k j t ) α ( 16
)
where A is a constant and α < 1.

2. A Cobb-Doublas specification with constant returns to scale:

f j (k j t , j t ) = A(k j t ) α ( j t ) 1-α ( 17
)
where α is the share of capital in the production function.

3. A constant elasticity of substitution (CES) function:

f j (k j t , j t ) = A α(k j t ) μ j + (1 -α)( j t ) μ j 1 μ j (18)
where 1 1-μ j is the elasticity of substitution between capital and labor.

Classification of variables and policy function

The model described above has the following characteristics: 12

• There are 4N + 1 endogenous (control) variables: consumption (c j t ), labor ( j t ), capital stock (k j t+1 ), investment (i j t ), Lagrange multiplier (λ t ).

• Consequently, there are 4N + 1 model equations given by ( 6)-(10).

• There are N exogenous variables (technology, a j t ), governed by equation (11).

• Among the endogenous variables, only previous period capital stock

k j
t is a predetermined state variable.

1 In models A1 and A5, there is no labor. There are therefore only 3N + 1 endogenous variables and as many model equations.

2 An alternative categorization of variables would include technology level a j t in the set of endogenous variables, thus making a total of 5N + 1 endogenous variables and as many model equations. The set of exogenous variables would then consist of the N + 1 shocks e j t and et, and there would be 2N predetermined endogenous state variables (k j t and a j t-1 ). We don't adopt this point of view since it artificially expands the dimension of the state space (it becomes of dimension 3N + 1).

• The state space is therefore of dimension 2N . 3 We denote by z t = (c t , t , i t , λ t ) the vector of the 3N + 1 endogenous variables (excluding capital stocks). The model can then be written in the following compact form:

E t H(z t+1 , z t , k t+1 , k t , a t+1 , a t ) = 0
The model solution is given by policy functions k t+1 = F (k t , a t ) and z t = G(k t , a t ), which satisfy:

E t H(G(F (k t , a t ), a t+1 ), G(k t , a t ), F (k t , a t ), k t , a t+1 , a t ) = 0

Specifications, values and deterministic steady state

We tried to select examples that are both representative of useful models and challenging from a numerical point of view. Finally, thirty different models were selected. Table 1 gives the values of the parameters which are constant across all specifications. Table 2 lists the various specifications with the functional form of the utility and production functions, and gives the values of the parameters which vary across specifications. 4 Note that, for parameter φ (the investment adjustment cost), we choose a low value, thereby increasing the variance of the endogenous variables: in test 2, which considers accuracy along a simulated time path, this leads to considering state points further away from the steady state.

The number of countries, N , considered in each model is different so as to check the ability of each algorithm to deal with both small and larger models. Some parameters still need to be specified: the Negishi weights τ j , and the parameters of the utility and production functions: A, b j (only for A2, A4, A6, A8), ψ j (only for A3, A7) and L e (only for A3, A4, A7, A8). These parameters are chosen so that: 3 As noted above, this relies on the fact that φ > 0. If there was no adjustment cost, the state space would be of dimension N + 1.

4 In this table, for asymmetric specifications A5-A8, the range notation Notes: N is the number of countries; u the equation number for the utility specification; f the equation number for the production function specification; γ, η, μ and χ, the value of the parameter when it appears in a given model specification.

Table 2: Model specifications

• At the steady state, all the countries consume exactly their net domestic production. This implies that in the dynamic model, net foreign assets remain small, which is a rough approximation of reality.

• At the steady state, all the countries have the same size. More precisely we assume that steady capital and labor supply are equal to unity (which implies that production of a country equals A, and therefore cj = A given the previous assumption).

• In all models, the steady state share of capital income in production is equal to α. This will be even true in A4 and A8, because we impose kj = ¯ j = 1.

• In models with a time endowment (A3, A4, A7, A8), the steady state ratio of labor supply over the time endowment is equal to 40%.

These assumptions lead to the following values for the remaining parameters:

A = 1 -β αβ τ j = 1 u j c ( cj , ¯ j ) = 1 u j c (A, 1) b j = (1 -α)A 1-1 γ j
(for A2, A6)

L e = 2.5 (for A3, A4, A7, A8)

b j = (1 -α) A 1-1 χ j (L e -¯ j ) -1 χ j = (1 -α) A 1-1 χ j (L e -1) -1 χ j
(for A4, A8)

ψ j = 1 L e -α(L e -¯ j ) = 1 L e -α(L e -1) (for A3, A7)
Given these restrictions on the parameters, the deterministic steady state for all the model specifications is:

¯ j = 1 (except A1, A5) kj = 1 āj = 1 īj = δ cj = A λ = 1
3 Accuracy tests

Approximation errors

We now turn to the definition of approximation errors. For given approximated policy functions ( F , Ĝ), we denote by R ( F , Ĝ) (k t , a t ) the vector of unit-free approximation errors at point (k t , a t ), determined as follows:

R j ( F , Ĝ) (k t , a t ) = τ j u j c (c j t , j t ) -λ t τ j u j c (c j t , j t ) R N +j ( F , Ĝ) (k t , a t ) = τ j u j (c j t , j t ) + λ t a j t f j (k j t , j t ) τ j u j (c j t , j t ) R 2N +j ( F , Ĝ) (k t , a t ) = λ t 1 + φ i j t k j t -δ -β E t λ t+1 1 + a j t+1 f j k (k j t+1 , j t+1 ) +φ 1 -δ + i j t+1 k j t+1 - 1 2 i j t+1 k j t+1 -δ i j t+1 k j t+1 -δ / λ t 1 + φ i j t k j t -δ R 3N +j ( F , Ĝ) (k t , a t ) = k j t+1 -(1 -δ)k j t -i j t k j t+1 R 4N +1 ( F , Ĝ) (k t , a t ) = N j=1 c j t + i j t -δk j t -a j t f j (k j t , j t ) + φ 2 k j t i j t k j t -δ 2 N j=1 c j t + i j t -δk j t -a j t f j (k j t , j t ) for j = 1, . . . , N, where k t+1 = F (k t , a t ), (c t , t , i t , λ t ) = G(k t , a t ) and (c t+1 , t+1 , i t+1 , λ t+1 ) = G(F (k t , a t ), a t+1 ).
By definition, the true policy function (F, G) is such that R (F,G) (k t , a t ) = 0 for all points of the state space.

Note that, from a computational point of view, evaluating R ( F , Ĝ) (k t , a t )

given F and Ĝ is straightforward, except for the approximation error of the Euler equation where an expectation term enters in; this requires the use of a numerical integration method, over the distribution of shocks which enter equation ( 11), as we describe it below.

Accuracy on a sphere (Test 1)

The first test compares the approximation errors of the participants' solutions on a hypersphere around the steady state of the model.

More formally, we denote x t = (k t , a t ) a point in the state space, x = ( k, ā) = (1, 1) the steady state and S r the hypersphere of radius r around the steady state given by S r = {x t ∈ R 2N ; x tx = r} (where is the Euclidian norm).

For the three values r = 0.01, r = 0.1 and r = 0.3, we draw 1,000 points in S r , and for every approximated policy function ( F , Ĝ), we compute the approximation errors R ( F , Ĝ) on these points, and report the maximum errors (for each equation separately, and for the aggregate across equations). This is the way errors are calculated in [START_REF] Judd | Projection methods for solving aggregate growth models[END_REF].

The purpose of this test is to quantify the loss of accuracy of a solution when going away from the steady state. It is expected that projection methods will perform better than perturbation methods for large values of the radius r, but this test will indicate to which extent this is true, and whether this is counterbalanced by a greater accuracy of perturbation methods at points very close to the steady state.

Accuracy along a simulated path (Test 2)

The second test compares the approximation errors along a simulated path, starting from the steady state of the model, and generating random draws of the shocks.

For a given approximated policy function ( F , Ĝ), we generate a series of points in the state space x 0 , . . . , x T where x t = (k t , a t ) and T =10,200.

The initial state is the steady state x 0 = ( k, ā), and the subsequent states are given by k t+1 = F (x t ) and the law of motion of a t (equation ( 11)).

The series of shocks (e 1 t , . . . , e N t , e t ) for t = 1, . . . , T are generated with a pseudo-random number generator.

The first 200 points in the generated series are discarded, so that the 10,000 remaining points can be considered as representative of the ergodic distribution of the model variables. On these points, we compute the approximation errors R ( F , Ĝ) , and report the mean and maximum errors (for each equation separately, and for the aggregate across equations). This is the way errors are calculated in [START_REF] Jin | Perturbation methods for general dynamic stochastic models[END_REF].

Den Haan-Marcet statistics (Test 3)

It was initially planned to use the Den Haan and Marcet (1994) statistics as a third testing device for participants' solutions.

The test statistics checks the orthogonality of the error residuals (as specified by R) with some arbitrary function of the state variables (called instruments), along a simulated path (like those used in accuracy test 2).

The instruments used here are a constant and the first and second order monomials of the state variables. For each approximated policy function, this exercise is repeated a given number of times, in order to obtain an empirical distribution of the test statistics. For the true policy function, the test statistics is distributed according to a χ 2 distribution. The accuracy criterion for an approximated policy function is therefore the closeness of the empirical distribution of the statistics to the χ 2 distribution.

It was finally decided to drop this test in the comparison project for the following two reasons:

• The computational time for the test is extremely large for high values of N (number of countries) and reasonable simulation lengths (T =10,000) when all the instruments are used. Computations are only feasible with a small number of countries, smaller simulation lengths or restricted sets of instruments.

• More importantly, the test seems unable to usefully discriminate between the participants' solutions. For the exercises which were performed, the distribution of the test statistics was very close across all solutions: they were either all very close to the χ 2 distribution, or all very far from it, depending on the simulation length and the set of instruments. Since this property was not extensively tested, it cannot parameterized functional form used for approximating the true solution.

By default the program computes test 1 and 2 for all the specifications (symmetric and asymmetric ones), and uses the same settings as those used in the comparison paper of [START_REF] Kollmann | Comparison of numerical solutions to a suite of multi-country real business cycle models[END_REF]. It is however possible to change the tasks to be performed (computing the Den Haan-Marcet statistics, restricting to only some solutions or model specifications), or to change some parameters of the tests (the number of points, the seed of the random number generators, the numerical integration method, or the normalization of the residual for the aggregate resource constraint

  The multi-country RBC models used in this comparison exercise are introduced in Den Haan et al. (2010). The four basic specifications try to introduce both a simple model such as the first one that considers fixed labor and a simple Cobb-Douglas function, and more complicated setups with endogenous labor, CES production functions and asymmetries between countries. The number of countries varies from two to ten. Altogether the comparison involves 30 different models.

  [a, b] for a parameter (say γ) means that γ j = a + (j-1)(b-a) N -1 . In other words, a (b) is the lowest (highest) value of γ across countries, and the N values for γ are evenly spread over [a, b].

Table 1 :

 1 Values of parameters common to all specifications

	Parameter Value
	β	0.99
	α	0.36
	δ	0.025
	σ	0.01
	ρ	0.95
	φ	0.5

be considered as a more general result, but nevertheless this apparent lack of discriminative power of the Den Haan-Marcet statistics for this exercise was deemed a sufficient reason to drop it.

Test bench

In order to test the participants' solutions on an equal basis, we wrote a testing program which performs the accuracy tests on all the solutions. With that test bench, one can be sure that:

• the random elements of the tests are the same across all solutions: for test 1, these random elements are the points on the sphere around the steady state; for test 2, they are the series of random shocks used to generate the simulated paths;

• when computing the approximation error of the Euler equation, the numerical integration method is the same across all solutions;

• the numerical precision used in the calculations is the same across all solutions.

Here we shortly describe some algorithmic choices made in the test bench, and then give an overview of the implementation.

Algorithmic choices

Three numerical integration methods are implemented in the test bench:

• the product four-point Gauss-Hermite quadrature, as described in Abramowitz and Stegun (1964, p. 890, eq. 25.4.46);

• the degree 5 monomial formula described in Judd (1998, p. 275, eq. 7.5.11);

• quasi-Monte Carlo integration, using Niederreiter quasi-random lowdiscrepancy sequence (see [START_REF] Bratley | Implementation and tests of low-discrepancy sequences[END_REF].

The default choice of the test bench is to use the Gauss-Hermite formula up to dimension 6 (i.e. for N ≤ 5 since there are N + 1 random shocks), and the monomial formula for higher dimensional problems (i.e. for N ≥ 6).

For test 1, the test bench generates uniformly distributed points on a hypersphere of dimension 2N using the simple algorithm described by [START_REF] Muller | A note on a method for generating points uniformly on n-dimensional spheres[END_REF]. The idea is to draw points from a multivariate standard Gaussian distribution of dimension 2N , and to divide these points by their Euclidian norm: it is easy to see that the resulting points are uniformly distributed over the unit hypersphere. Points uniformly distributed over any hypersphere can immediately be generated by rescaling and translating. The Gaussian draws are generated using draws from the Sobol quasi-random low-discrepancy sequences (see [START_REF] Antonov | An economic method of computing LP[tau]-sequences[END_REF] which are transformed in Gaussian vectors by applying the reciprocal of the Gaussian cumulative distribution function.

For test 2, the random draws for shocks are generated using a Mersenne-Twister random number generator (see [START_REF] Matsumoto | Mersenne twister: a 623dimensionally equidistributed uniform pseudo-random number generator[END_REF].

Implementation choices

The test bench is a standalone program written in the C++ programming language, and relies on the GNU Scientific library 5 for linear algebra computations and random number generation. The source code of the program is freely available on the web. 6

The program is made of two main modules, plus one module per solution to be tested. The first module implements the 30 specifications and computes the residuals. 7 The second one implements the accuracy tests. 8

Then for each of the six participants' solutions (PER1, PER2, MRGAL, SMOL, CGA and SSA1, using the notations of [START_REF] Kollmann | Comparison of numerical solutions to a suite of multi-country real business cycle models[END_REF], a module implements the corresponding approximated policy function. 9 Note that the test bench does not contain the code to fully recompute each solution: 10 it is only able to simulate the solution, using the coefficients of the 5 See Galassi et al. (2003) or http://www.gnu.org/software/gsl. 6 The source code that was used for generating the comparison results presented in [START_REF] Kollmann | Comparison of numerical solutions to a suite of multi-country real business cycle models[END_REF] is available as supplemental material to this paper, and is downloadable at: .... Future versions of the software will be available at: http: //www.dynare.org/JedcTestsuiteWiki/TestBenchProblemA.

7 See class ModelSpec and its subclasses ModelSpecA1 to ModelSpecA8.

8 See class SolutionTester and main file tester.cc. 9 Note that SMOL, CGA and SSA1 approximated policy functions do not return a value for the Lagrange multiplier λt. The workaround is the following: for SMOL solution, the value of τ 1 u 1 c (c 1 t , 1 t ) is used as a replacement for λt; for CGA and SSA, the arithmetic mean of τ j u j c (c j t , j t ) over j = 1, . . . , N is used as a replacement for λt. 10 The source code for recomputing the various solutions is available as supplemental


Roland Tóth
email: rolandtoth@berkeley.edu

Vincent Laurain
email: vincent.laurain@cran.uhp-nancy.fr

Marion Gilson
email: marion.gilson@cran.uhp-nancy.fr

Hugues Garnier
email: hugues.garnier@cran.uhp-nancy.fr

Instrumental variable scheme for closed-loop LPV model identification

Keywords: Linear parameter-varying system, Identification, Refined instrumental variable, Box-Jenkins models

Identification of real-world systems is often applied in closed loop due to stability, performance or safety constraints. However, in case of Linear Parameter-Varying (LPV) systems, closed-loop identification is not well-established despite the recent advances in prediction error approaches. Building on the available results, the paper proposes the closed-loop generalization of a recently introduced instrumental variable scheme for the identification of LPV-IO models with Box-Jenkins type of noise model structures. Estimation under closed-loop conditions with the proposed approach is analyzed from the stochastic point of view and the performance of the method is demonstrated through a representative simulation example.

Introduction

Identification of physical or chemical systems is often restricted to data gathered during closed-loop operation due to stability, performance/economical or safety constraints. As the control loop itself introduces correlation between the disturbances and the control signal, the statistically optimal (unbiased with minimal variance) estimation of the parameters of a chosen model structure w.r.t. the data-generating system is an essentially different problem than in the open-loop setting. Hence in the identification literature, especially in the linear timeinvariant (LTI) case, many approaches have been proposed to provide well-applicable solutions in this problem setting. An overview of the recent developments in the LTI case can be found in [START_REF] Van Den Hof | Closed-loop issues in system identification[END_REF] and [START_REF] Forssell | Closed-loop identification revisited[END_REF].

Identification of linear parameter-varying (LPV) models for systems operated in open loop has recently seen significant improvements. The case of closed-loop LPV model identification has however remained sparsely studied (see [START_REF] Casella | LPV/LFT modelling and identification: overview, synergies and a case study[END_REF][START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]). The main difficulty in this system class in terms of identification is that even if the signal relations between the inputs u and outputs y of the system are linear, these relations are allowed to change over time as a function of a so-called scheduling variable p : Z → P with P ⊆ R np . This allows to describe a large class of nonlinear/time-varying systems in an attractive structure on which the well-established LPV controlsynthesis approaches, e.g. [START_REF] Packard | Gain scheduling via linear fractional transformations[END_REF][START_REF] Scherer | Mixed H 2 /H∞ control for time-varying and linear parametrically-varying systems[END_REF], can be applied. On the other hand, this varying-relation prevents the use of crucial system theoretical relations, like transfer functions and commutativity of operators, on which most LTI closed-loop identification approaches are based. As a consequence, only preliminary closed-loop methods have been proposed in the literature without being able to exploit the existing tools and knowledge available in the LTI case. In [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF], an approximation based LPV extension of a predictor subspace approach (PSBID) has been introduced which is also applicable in a closedloop setting. In [START_REF] Boonto | Closed-loop system identification of LPV input-output models -application to an arm-driven pendulum[END_REF], also an approximation based LPV extension of the CLOE algorithm (see [START_REF] Landau | Recursive algorithms for identification in closed loop: A unified approach and evaluation[END_REF]) has been investigated w.r.t. LPV output-error (OE) type of models. In [START_REF] Abbas | An instrumental variable technique for open-loop and closed-loop identification of input-output LPV models[END_REF] and [START_REF] Butcher | On the consistency of certain identification methods for linear parameter varying systems[END_REF], a basic instrumental variable method has been introduced. It is a clear indicator of the immature state of this research direction that the stochastic properties of the estimation has only been analyzed in [START_REF] Butcher | On the consistency of certain identification methods for linear parameter varying systems[END_REF]. Unfortunately in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF], it has been shown recently that the formulation of the one-step IV approach, proposed in [START_REF] Butcher | On the consistency of certain identification methods for linear parameter varying systems[END_REF], does not allow to reach statistically optimal esti-mates due to the non-commutativity of certain filtering operations. This highlights that up to now, no method has been established which allows statistically optimal estimation under a closed-loop setting especially with general noise models like Box-Jenkins (BJ).

Based upon the available results of [START_REF] Abbas | An instrumental variable technique for open-loop and closed-loop identification of input-output LPV models[END_REF][START_REF] Butcher | On the consistency of certain identification methods for linear parameter varying systems[END_REF] and using the theoretical advancements of the LPV prediction errorframework introduced in [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF], we propose in this paper a closed-loop extension of a recently developed IV approach for LPV-BJ models [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF]. In this setting, the noise v o affecting the sampled output measurement y(t k ) is assumed to have a rational spectral density which is allowed to be unrelated to the actual process dynamics of the data-generating system (general noise structure). As a first step towards the case of a p-dependent noise, it is also assumed that this rational spectral density does not depend on p. A stochastic analysis of the proposed closed-loop approach is provided, exploring the limitations and the benefits of this estimation strategy. The performance of the algorithm is also demonstrated on a representative simulation example.

Problem description

Consider the data-generating LPV system S o given in the closed-loop setting of Fig. 1, where S o is defined as:

A o (p k , q -1 )χ o (t k ) = B o (p k , q -1 )u(t k-d ), (1a) y(t k ) = χ o (t k ) + v o (t k ). (1b)
Here u(t k ) is the input of the plant, p k is the value of the scheduling variable p at sample time t k , χ o is the noise-free output, v o is the additive noise with bounded spectral density, y is the noisy output of the system, d is the delay, and q is the time-shift operator, i.e. q -i u(t k ) = u(t k-i ). A o (p k , q -1 ) and B o (p k , q -1 ) are polynomials in q -1 of degree n a and n b respectively:

A o (p k , q -1 ) = 1 + na i=1 a o i (p k )q -i , (2a)
B o (p k , q -1 ) = n b j=0 b o j (p k )q -j , (2b)
where the coefficients a o i and b o j are real meromorphic functions2 with static dependence on p, i.e. dependence only on the instantaneous value of p at time t k . It is assumed that each a o i and b o j is non-singular on P, thus the solutions of S o are well-defined and the process part G o , defined by (1a), is completely characterized by the coefficient functions {a o i } na i=1 and {b o j } n b j=0 . The noise v o is assumed to be independent from p. Latter we will return to the more general case of p-dependent noise models, showing that in case of an IV-based estimation, violation of this assumption does not lead to a biased estimate.

r 2 (t k ) G o C o r 1 (t k ) u(t k ) p k e o (t k ) v o (t k ) χ o (t k ) y(t k ) H o (q -1 )
Fig. 1. Closed-loop LPV system.

If v o has a rational spectral density, then the noise process can be represented by an autoregressive moving average (ARMA) model:

v o (t k ) = H o (q)e o (t k ) = C o (q -1 ) D o (q -1 ) e o (t k ), (3)
where C o (q -1 ) and D o (q -1 ) = 0 are monic polynomials with constant coefficients and with respective degree n c and n d . The corresponding proper transfer function H o (q) is assumed to be stable and to have a stable inverse. In case C o (q -1 ) = D o (q -1 ) = 1, (3) defines an OE noise model, however with C o (q -1 ) = D o (q -1 ), ( 3) is general enough to represent BJ-type of noise models.

In terms of the closed-loop setting of Fig. 1,

u(t k ) = r 1 (t k ) + C o (r 2 (t k ) -y(t k ))
, where C o is the operator form of the controller and r 1 , r 2 are reference signals. Note that, as in many digital control systems, u is known (the error introduced by the actuators can be modeled as an output additive (p-dependent) noise due to the linearity of the system). The controller can be any LTI, nonlinear or LPV controller, under the assumptions that • C o is a priori known;

• C o ensures BIBO stability of the closed-loop system S o for any p ∈ P Z .

Next we introduce a model structure and parametrization for the identification of S o , where, according to (1ab) and ( 3), the noise model and the process model are parameterized separately. The proposed LPV-BJ model, denoted in the sequel as M θ , is defined as:

A(p k , q -1 , ρ)χ(t k ) = B(p k , q -1 , ρ)u(t k-d ), (4a) D(q -1 , η)v(t k ) = C(q -1 , η)e(t k ), (4b) y(t k ) = χ(t k ) + v(t k ), ( 4c
) with parameters θ = [ ρ ⊤ η ⊤ ] and with u(t k ) = r 1 (t k )+ C o (r 2 (t k ) -y(t k )).
The process model part of M θ , denoted by G ρ , is defined in terms of the LPV-IO representation (4a) where A and B are polynomials with order n a and n b respectively and with p-dependent coefficients a i and b j parameterized as

a i (p k ) = a i,0 + nα l=1 a i,l f l (p k ), i = 1, . . . , n a , ( 5a
) b j (p k ) = b j,0 + n β l=1 b j,l g l (p k ), j = 0, . . . , n b . (5b)
In this parametrization, {f l } nα l=1 and {g l } n β l=1 are a priori chosen meromorphic functions of p bounded on P, with static dependence, allowing the identifiability of the model (e.g. linearly independent functions on P). The associated model parameters ρ are stacked columnwise:

ρ = [ a 1 . . . a na b 0 . . . b n b ] ⊤ ∈ R nρ where a i = [ a i,0 a i,1 . . . a i,nα ], b j = [ b j,0 b j,1 . . . b j,n β ] and n ρ = n a (n α + 1) + (n b + 1)(n β + 1).
The noise-model part of M θ , denoted by H η , is defined in terms of (4b) where C and D are monic polynomials with order n c and n d respectively and with constant coefficients c i and

d j collected as η = [ c 1 . . . c nc d 1 . . . d n d ] ⊤ ∈ R nη , where n η = n c + n d . As (4b) is LTI, it can be represented by a transfer function H(q, η) = C(q -1 , η)/D(q -1 , η). It is further assumed that η ∈ Θ η ⊂ R nη such that H(z, η)
is stable (analytic in the exterior of the unit circle) and has a stable inverse.

Introduce also G = {G ρ | ρ ∈ R nρ } and H = {H η | η ∈ R nη },
as the collection of all process and noise models in the form of (4a) and (4b). Then the model set, denoted as M, takes the form

M = (G ρ , H η ) | col(ρ, η) = θ ∈ R nρ+nη , (6)
corresponding to the set of candidate models in which we seek the model that explains data gathered from S o the best under a given identification criterion.

Denote D N = {y(t k ), u(t k ), p(t k )} N k=1 a data sequence of S o . In [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF], it has been shown, that a one-step-ahead output predictor ŷθ (t k |t k-1 ) can be formulated w.r.t. the considered model structure (4a-c) under the commonly used assumption that noise-free observation of the sequence {p k , p k-1 , . . .} is available. Recently it has been proven that using estimated moments, such a predictor can be formulated if p is observed up to an additive white noise independent from v o [START_REF] Tóth | A prediction-error identification framework for linear parameter-varying systems[END_REF]. However for the sake of simplicity, we will only consider the previous case. These results allow to formulate (in either case) the estimation of θ, as the minimization of the one-step ahead prediction error :

ε θ = y(t k ) -ŷθ (t k |t k-1
). This minimization is formulated in terms of an identification criterion W (D N , θ), like the least squares criterion

W (D N , θ) = 1 N N k=1 ε 2 θ (t k ) = 1 N ε θ 2 ℓ2 , (7)
such that the parameter estimate is θN = arg min

θ∈R nρ+nη W (D N , θ). (8)
Based on the previous considerations, the identification problem addressed in the sequel is defined as follows:

Problem 1 Given a discrete-time LPV data-generating system S o in the closed-loop setting of Fig. 1 with a priori known stabilizing controller C o . Based on the LPV-BJ model structure M θ defined by (4a-c) and a data set D N collected from S o , estimate the parameter vector θ as the minimization of (7) under the following assumptions:

A1 S o ∈ M. A2 In the parametrization (5a-b) of A and B, {f l } nα l=1 and {g l } n β l=1 are chosen such that M θ is identifiable. A3 D N is informative w.r.t. M θ : if θ 1 = θ 2 then M θ1
and M θ2 lead to different prediction errors given D N . A4 S o is BIBO stable, i.e. for any bounded p ∈ P Z and u ∈ R Z , the output of S o is bounded [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]. A5 The closed-loop system with C o is BIBO stable.

In [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF], a refined instrumental variable (RIV) based approach has been introduced to solve the estimation problem of M θ in an open-loop setting. However, a major difference between the open-loop setting considered in [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF] and the closed-loop setting of Prob. 1 is the correlation of u with e o due to the feedback loop. Therefore in the following, we aim to solve the problem of determining a suitable instrument that allows to developing the LPV-RIV approach in the closed-loop case.

RIV approach for closed-loop LPV systems

In [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF] it has been shown that reformulation of (4a-c) as a multiple-input single-output (MISO) LTI model allows an elegant solution of the filtering problem associated with LPV-IV approaches. This reformulation is necessary to avoid the problem of non-commutativity of parameter-varying filters and hence required to solve Prob. 1. Consequently, (4a) is rewritten as

χ(t k ) + na i=1 a i,0 χ(t k-i ) F (q -1 )χ(t k ) + na i=1 nα l=1 a i,l f l (p k )χ(t k-i ) χ i,l (t k ) = n b j=0 n β l=0 b j,l g l (p k )u(t k-d-j ) u j,l (t k ) , (9)
where F (q -1 ) = 1 + na i=1 a i,0 q -i and g 0 ( ) ≡ 1. Note that in this way, the LPV-BJ model is rewritten as a MISO system with (n b + 1)(n β + 1) + n a n α inputs {χ i,l } na,nα i=1,l=1 and {u j,l } n b ,n β j=0,l=0 . F (q -1 ) does not depend on p k , thus ( 9) and (4b-c) have the following LTI form:

y(t k ) = - na i=1 nα l=1 a i,l F (q -1 ) χ i,l (t k ) + n b j=0 n β l=0 b j,l F (q -1 )
u k,j (t k ) + H(q)e(t k ). [START_REF] Ljung | System Identification, theory for the user[END_REF] Note that ( 10) is equivalent with (4a-c), but it is not a representation of the associated LPV system as it includes the lumped output variables in {χ i,l } na,nα i=1,l=1 . Using [START_REF] Ljung | System Identification, theory for the user[END_REF], the estimation problem of the parameters can be formulated as a linear regression which allows optimal solution of (7) w.r.t. LPV-BJ models. To achieve this solution, one possible way is to develop an extension of the closed-loop refined instrumental variable (RIV) approach of the LTI framework. Next we derive this extension, which, as we will see, provides an easily implementable iterative estimation scheme.

Rewrite [START_REF] Ljung | System Identification, theory for the user[END_REF] to the linear regression form:

y(t k ) = ϕ ⊤ (t k )ρ + ṽ(t k ), (11)
where ṽ(t k ) = F (q -1 , ρ)v(t k ) and

ϕ(t k ) = [ -y(t k-1 ) . . . -y(t k-na ) -χ 1,1 (t k ) . . . -χ na,nα (t k ) u 0,0 (t k ) . . . u n b ,n β (t k ) ] ⊤ , ρ = [ a 1,0 . . . a na,0 a 1,1 . . . a na,nα b 0,0 . . . b n b ,n β ] ⊤ .
Two difficulties still remain to obtain the minimum of (7) based on [START_REF] Packard | Gain scheduling via linear fractional transformations[END_REF]: the regressor ϕ(t k ) contains the unknown terms {χ i,l (t k )} na,nα i=1,l=0 and all of its elements are corrupted with the colored noise v(t k ). To resolve this problem, an appropriate instrument ζ(t k ) can be introduced such that the estimate of ρ can be given as [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]:

ρXIV (N ) = arg min ρ∈R nρ 1 N N k=1 L(q)ζ(t k )L(q)ϕ ⊤ (t k ) ρ - 1 N N k=1 L(q)ζ(t k )L(q)y(t k ) 2 W ( 12
)
where x 2 W = x T W x, with W a positive definite weighting matrix and L(q) is a stable pre-filter. If G o ∈ G, the estimate ( 12) is consistent under the following well-known conditions 3 :

C1 Ē{L(q)ζ(t k )L(q)ϕ ⊤ (t k )} is full column rank. C2 Ē{L(q)ζ(t k )L(q)ṽ(t k )} = 0.
Moreover, the minimum variance estimator can be achieved if [START_REF] Gilson | Instrumental variable methods for closed-loop system identification[END_REF][START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF][START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF]:

C3 W = I and ζ(t k
) is chosen as the noise-free version of [START_REF] Packard | Gain scheduling via linear fractional transformations[END_REF] and is therefore defined in the LPV case as:

ζ(t k ) = -χ(t k-1 ) . . . -χ(t k-na ) -χ 1,1 (t k ) . . . -χ na,nα (t k ) ů0,0 (t k ) . . . ůn b ,n β (t k ) ⊤
where ů and χ are the signals from the auxiliary model as presented in Fig. 2. C4 G o ∈ G and n ρ is equal to the minimal number of parameters required to represent G o with the considered model structure. 3 The notation Ē{.} = limN→∞ 1 N N t=1 E{.} is adopted from the prediction error framework of [START_REF] Ljung | System Identification, theory for the user[END_REF].

r 2 (t k ) C o r 1 (t k ) ů(t k ) p k χ(t k ) G o
Fig. 2. Auxiliary model C5 L(q) is chosen as [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF]:

L(q) = D(q -1 ) F (q -1 )C(q -1 ) . ( 13
)
Full column rank of Ē{L(q)ζ(t k )L(q)ϕ ⊤ (t k )} is ensured by the informativity of the data set (A3). Note that conditions about informativity w.r.t. LPV-BJ model structures is in the focus of current research (see [START_REF] Bamieh | Identification of linear parameter varying models[END_REF][START_REF] Butcher | On the consistency of certain identification methods for linear parameter varying systems[END_REF][START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF] for results in the ARX case). To fulfill C1 under A3, the discussion can be found in [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]. It is also important to note that in a practical situation none of F (q -1 , ρ), C(q -1 , η), D(q -1 , η) or {a i,l (ρ)} na,nα i=1,l=0 , {b j,l (ρ)} n b ,n β j=0,l=0 is known a priori. Therefore, the RIV estimation normally involves an iterative (or relaxation) algorithm in which, at each iteration, an 'auxiliary model' is used to generate the instrumental variables (which guarantees C2), as well as the associated pre-filters. Convergence of such an iterative RIV algorithm has not been proven so far and is only empirically assumed [START_REF] Young | The refined instrumental variable method: Unified estimation of discrete and continuous-time transfer function models[END_REF].

Based on the previous considerations, the iterative scheme of the RIV algorithm can be extended to the closed-loop LPV case as described by Algorithm 1.

Based on a similar concept, the so-called simplified LPV-CLRIV (LPV-CLSRIV) method, can also be developed for the estimation of LPV-OE models. This method is based on a model structure (4a-c) with C(q -1 , η) = D(q -1 , η) = 1 and consequently, Step 7 of Algorithm 1 can be skipped.

Remark 1 During the formulation of the approach it has been assumed that v o is independent from p (see (3)). However, to have a general noise model in the LPV setting it is necessary to consider p-dependence of the noise as has been pointed out in [START_REF] Tóth | Modeling and Identification of Linear Parameter-Varying Systems[END_REF]. However the RIV setting, even using an OE type of noise model (like in the SRIV case), leads to unbiased estimates also in the case when the assumed noise model is incorrect (or correlated with p). The role of a more accurate noise model in this setting is to lower the variance of the resulting estimates. This explains why using a general BJ structure is important to increase the accuracy of the RIV method. However in the LPV case, allowing the noise model to be also p-dependent not only makes the estimation problem more difficult from the filtering point of view, but also increases significantly the number of parameters to be estimated. This results in a trade-off problem, where the simple estimation scheme proposed in this paper provides an attractive and practically applicable approach. Finding an Algorithm 1 LPV-CLRIV 1: set τ = 0 and let θ(0) = [ (ρ (0) ) ⊤ (η (0) ) ⊤ ] ⊤ be an initial parameter estimate given by the least square (LS) approach (ARX estimate of M θ ). 2: repeat 3:

compute an estimate of χ(t k ) via

A(p k , q -1 , ρ(τ) ) χ(t k ) = B(p k , q -1 , ρ(τ) )u(t k-d ),
where G ρ(τ) is assumed to be stable (A5) and compute { χi,l (t k )} na,n α i=1,l=0 according to [START_REF] Laurain | Refined instrumental variable methods for identification of LPV Box-Jenkins models[END_REF].

4:

compute the estimated filter:

L(q -1 , θ(τ) ) = D(q -1 , η(τ) ) C(q -1 , η(τ) )F (q -1 , ρ(τ) )
and the filtered signals

{u f j,l (t k )} n b ,n β j=0,l=0 , y f (t k ) and { χf i,l (t k )} na,nα i=1,l=0 . 5:
build the estimated filtered regressor as:

φf (t k ) = [ -y f (t k-1 ) . . . -y f (t k-na ) -χf 1,1 (t k ) . . . -χf na,nα (t k ) u f 0,0 (t k ) . . . u f n b ,n β (t k ) ] ⊤
and compute the filtered instrument ζf (t k ) by simulating G ρ(τ) according to Fig. 2:

ζf (t k ) = [ -χf (t k-1 ) . . . -χf (t k-na ) -χf 1,1 (t k ) . . . -χf na,nα (t k ) ûf 0,0 (t k ) . . . ûf n b ,n β (t k ) ] ⊤ 6:
With W = I and L(q -1 , θ(τ) ) (according to C3-C5), compute the solution of (12) via

ρ(τ+1) = 1 N N k=1 ζf (t k ) φ⊤ f (t k ) -1 N k=1 1 N ζf (t k )y f (t k ). 7:
compute the estimate of the noise signal v as

v(t k ) = y(t k ) -χ(t k , ρ(τ) ).
Based on v, estimation of η(τ+1) follows using e.g. the ARMA estimation algorithm of the MATLAB identification toolbox (an IV approach can also be used for this purpose, see [START_REF] Young | The refined instrumental variable method: Unified estimation of discrete and continuous-time transfer function models[END_REF]).

8:

increase τ by 1 9: until θ (τ ) has converged IV-based estimation scheme for p-dependent BJ models, which preserves the simplicity of the proposed approach, could lead to further decrease of the variance of the estimates and hence it is in the focus of future research.

Simulation Example

In the following example we aim to demonstrating the performance of proposed closed-loop RIV scheme on a relevant simulation example. In this demonstration, we will focus on estimation under different noise conditions/models, hence the model structure of the process part is assumed to be known.

Consider the example of a mass connected to a spring and a varying damper depicted in Fig. 3. This problem is one of the typical phenomena occurring in the motion control of many mechatronic systems like in active suspension. Denote x the position (in [m]) of the mass m (in [kg]), k s > 0 the stiffness of the spring and c d > 0 the varying damping. Introduce F as the force (in [N]) acting on the mass m. Then in continuous time (CT), the behavior of the system is defined by:

d dt m d dt x(t) = F (t) -k s x(t) -c d (t) d dt x(t). (14)
By considering F (t) as the input u(t), x(t) as the output y(t) and p(t) : R → [0, 1] as the scheduling variable such that c

d (t) = c (0) d +c (1)
d p(t), ( 14) can be rewritten as a CT-LPV system. By using a simple backward Euler type of discretization in a zero-order-hold setting with sampling period T d > 0, this LPV system can be formulated as a discrete-time (DT) LPV system:

y(t k ) = T d c (1) d p k +T d c (0) d +2m L(p k ) y(t k-1 ) -m L(p k ) y(t k-2 ) + T 2 d L(p k ) u(t k ) (15)
where u(t k ), y(t k ), p k denote the sampled signals of F (t), x(t), and p(t) respectively and

L(p k ) = T d c (1)
d p k + T d c (0) d + m + k s T 2 d .
To simplify the problem, approximate 1 L by its 1 st -order Taylor approximation at the mid-point of P, i.e. at p k = 0.5:

1 L(p) p=0.5 ≈ 1 τ1 2 + τ 0 - τ 1 ( τ1 2 + τ 0 ) 2 (p -0.5), (16)
where

τ 0 = T d c (0) d + m + k s T 2 d and τ 1 = T d c (1)
d . Then the resulting DT-LPV representation reads as [START_REF] Van Wingerden | Subspace identification of bilinear and LPV systems for open-and closed-loop data[END_REF] where

y(t k ) = -a 10 + a 11 p k + a 12 p 2 k y(t k-1 ) - a 20 + a 21 p k y(t k-2 ) + b 00 + b 01 p k u(t k ),
a 11 = -τ ′ T d c (1) d + τ ′′ (T d c (0) d + 2m), a 10 = -τ ′ (T d c (0) d + 2m), a 12 = τ ′′ T d c (1) d , a 20 = τ ′ m, a 21 = -τ ′′ m, b 00 = τ ′ T 2 d , b 01 = -τ ′′ T 2 d with τ ′ = 2 τ1+2τ0 + 2τ1 (τ1+2τ0) 2 and τ ′′ = 4τ1 (τ1+2τ0) 2 .
Consider [START_REF] Tóth | A prediction-error identification framework for linear parameter-varying systems[END_REF] as G o in the closed loop setting of Fig. 1 with parameters T d = 0.05s, m = 0.01, k s = 0.85, c

= 0.5, c (1) d
= 0.5. To control the motion of this system an LPV-PI controller C o has been designed:

u(t k ) = u(t k-1 ) + K C (p k ) w(t k ) -w(t k-1 ) + K C (p k )T d T I (p k ) w(t k ), (18)
S o        A o (q, p k ) = 1 + a o 1 (p k )q -1 + a o 2 (p k )q -2 B o (q, p k ) = b o 0 (p k ) H o (q) = 1 1 -0.5q -1 -0.2q -2 (19) where v(t k ) = H o (q)e o (t k ), e o (t k ) ∈ N (0, σ 2 e ) and a o 1 (p k ) = -1.135 -0.1738p k + 0.2538p 2 k , (20a) a o 2 (p k ) = 0.2523 -0.1015p k , (20b) b o 0 (p k ) = 0.0631 -0.0254p k . (20c)
In the sequel, both the scheduling signal p and the input r(t k ) are considered as white noises with uniform distributions U(0.05, 0.95) and U(-0.5, 0.5) respectively. Realizations of these signals with length N = 2000 are used to generate data sets D N of S o . To demonstrate the efficiency of the presented method, the proposed approach is compared to the conventional least-squares (LPV-LS) method (see [START_REF] Bamieh | Identification of linear parameter varying models[END_REF]) and also to the direct nonlinear minimization of the least-squares criteria (7) (NL-PEM) w.r.t. the LPV-BJ model. The latter approach is computed via the LSQNONLIN algorithm of the optimization toolbox in MATLAB. Note that the NL-PEM approach provides statistically optimal results in case the model structure is optimally chosen and the nonlinear optimization succeeds to find the global minimum of (7) (depending on the initialization). Therefore to use the NL-PEM as a base-line result for assessing the performance of the IV approach, it is initialized with the true parameters and hence the provided estimates can be considered optimal.

The LPV-LS method assumes an LPV-ARX type of structure defined in this case as:

M LPV-ARX θ    A(p k , q -1 , ρ) = 1 + a 1 (p k )q -1 + a 2 (p k )q -2 B(p k , q -1 , ρ) = b 0 (p k ) H(p k , q, ρ) = A † (p k , q -1 , ρ)
where A † is the adjoint of A and

a 1 (p k ) = a 1,0 + a 1,1 p k + a 1,2 p 2 k , (21a) a 2 (p k ) = a 2,0 + a 2,1 p k , (21b) b 0 (p k ) = b 0,0 + b 0,1 p k . (21c)
The proposed closed-loop LPV Refined Instrumental Variable method (LPV-CLRIV) and the NL-PEM method are based on the LPV-BJ model structure: Table 1 shows the norm of the bias (BN) ||ρ o -Ē(ρ)|| ℓ2 and variance norm (VN) || Ē(ρ -Ē(ρ))|| ℓ2 of the estimated parameter vector w.r.t. the process part, where Ē is the mean operator over the Monte Carlo simulation. The table also displays the mean number of iterations (Nit) the algorithms needed to converge to the estimated parameter vector. Table 1 demonstrates that the LPV-CLRIV method is unbiased according to the theoretical results, while the LPV-LS method exhibits a considerable bias. For SNR down to 8dB, the LPV-CLRIV achieves a lower variance for ρ than the LPV-LS method. Moreover, the achieved variance is close to result of the statistically optimal NL-PEM estimates. To show the advantage of the proposed IV method w.r.t. NL-PEM approach, consider a more realistic scenario when the model structure of the noise is unknown. Thus, identification of S o is considered with the following LPV-OE model structure (corresponding to G o ∈ G, H o / ∈ H):

M LPV-BJ θ        A(p k , q -1 , ρ) = 1 + a 1 (p k )q -1 + a 2 (p k )q -2 B(p k , q -1 , ρ) = b 0 (p k ) H(p k , q, η) = 1 1 + d 1 q -1 + d 2 q -2 with a 1 (p k ), a 2 (p k ), b 0 (p k ), b 1 (p k )
M LPV-OE θ   
A(p k , q -1 , ρ) = 1 + a 1 (p k )q -1 + a 2 (p k )q -2 B(p k , q -1 , ρ) = b 0 (p k )q -1 + b 1 (p k )q -2 H(p k , q, η) = 1

This model is identified using the proposed simplified closed-loop RIV approach (LPV-CLSRIV) and the NL-PEM method. The results are presented in Table 2. It appears that the NL-PEM is affected by the incorrect structure of the noise model both in terms of bias and variance, while the proposed LPV-CLSRIV method only suffers from a mild variance increase (remaining close to the LPV-LS variance) but remains unbiased. Consequently, this representative example suggests that the proposed method is robust to the modeling error on the noise and displays a better bias/variance trade-off.

Conclusion

In this paper, a closed-loop LPV-RIV approach is introduced to provide an efficient solution for the closed-loop identification of LPV systems with Box-Jenkins type of noise models. The approach is formulated with the assumption of prior knowledge of the controller but without any restriction on its structure. It is shown that under given conditions the proposed method provides consistent estimates. The performance of the approach is demonstrated on a representative example pointing out that its particular advantage is being robust to noise and modeling errors.

Fig. 3 .

 3 Fig. 3. Mass-spring-damper system with varying damping.wherew(t k ) = r(t k )y(t k ), K C (p k ) = 3.33 + 0.99p k + 2.06p 2k and T I (p k ) = 0.225 + 0.215p k -0.09p 2 k . Additionally, a BJ type of noise v is added to the output. Consequently the data generating system can be expressed:

  as given in (21a-c), and hence represent the situation S o ∈ M. The robustness of the proposed and existing algorithms are investigated with respect to different signal-to-noise ratios: SNR = 10 log P χo Pv o , where P χo and P vo are the average power of the signals χ o and v o respectively. To provide representative results, a Monte Carlo simulation of N MC = 200 runs with new noise realizations is accomplished at SNR levels: 18dB, 13dB and 8dB.

Table 1

 1 Bias and variance results using the full model structure M LPV-BJ

	θ	.
		Method		18dB	13dB	8dB
		LPV-LS	BN	0.0849 0.2158 0.4134
			VN	0.0752 0.1232 0.1876
		NL-PEM	BN	40e-04 71e-04 0.0138
			VN	0.0464 0.0819 0.1393
		LPV	BN	39e-04 62e-04 0.0103
		CLRIV	VN	0.0505 0.0922 0.1638
			Nit	9	12	12

Table 2

 2 Bias and variance results using the inaccurate model structure M LPV-OE

	θ	.
	Method		18dB	13dB	8dB
	LPV-LS	BN	0.0849 0.2158 0.4134
			VN	0.0752 0.1232 0.1876
	NL-PEM	BN	0.010	0.027	0.072
			VN	0.0621 0.1069 0.1742
		LPV	BN	34e-04 59e-04 0.0134
	CLSRIV	VN	0.0612 0.1085 0.2124
			Nit	13	15	18

f : R n → R is a real meromorphic function if f = g/h with g, h analytic and h = 0.

(Hugues Garnier). 1 Supported by the Netherlands Organization for Scientific Research (grant no. 680-50-0927).


Jean-Marc Chassery
email: jean-marc.chassery@gipsa-lab.grenoble-inp.fr

Isabelle Sivignon
email: isabelle.sivignon@gipsa-lab.grenoble-inp.fr

Optimal covering of a straight line applied to discrete convexity

The relation between a straight line and its digitization as a digital straight line is often expressed using a notion of proximity. In this contribution, we consider the covering of the straight line by a set of balls centered on the digital straight line pixels. We prove that the optimal radius of the balls is strictly less than one, and can be expressed as a function of the slope of the straight line. This property is used to define discrete convexity in concordance with previous works on convexity.

Introduction

From the seminal work of Sklansky [START_REF] Sklansky | Recognition of convex blobs[END_REF], discrete convexity has been the subject of many studies with a common objective: transcribe the Euclidean definition in the digital space.

In the eighties, a previous work introduced ε-convexity using covering of connected sets by balls [START_REF] Chassery | Discrete convexity: Definition, parametrization, and compatibility with continuous convexity[END_REF]. In the case of convex shapes, it was shown in [START_REF] Chassery | Discrete convexity: Definition, parametrization, and compatibility with continuous convexity[END_REF] that the value ε can be written as p p+1 (lower than 1) where p is a parameter computed from the edges of the convex hull of the shape. However, the value of this parameter p was given algorithmically but not analytically. In this paper, using the arithmetical definition of a digital straight line [START_REF] Reveilles | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF][START_REF] Debled-Rennesson | A linear algorithm for segmentation of digital curves[END_REF], we prove that ε can be expressed exactly as a function of the characteristics of the digital straight line supporting the edges of the convex hull of the shape.

After classical definitions of digital straight line and its characteristics, we establish the property of covering a straight line by balls centered exclusively on digital straight line pixels. The following section applies that property to the case of discrete convexity and we present the algorithm for discrete convex hull computation working only on digital space.

Preliminary definitions

Let L be an Euclidian straight line in R 2 given by the equation axby + µ = 0, with a, b, µ in Z, and gcd(a, b) = 1. In the following, we also assume without loss of generality that 0 ≤ a ≤ b. All other cases are symmetrical. Such a straight line may be considered as supporting a linear contour of a shape, such that all the points of the shape are on the same side of the straight line.

Let us now consider the Object Boundary Quantization of L on the isothetic regular grid. It is given by the set of pixels (x, y) such that x ∈ Z and y = -ax-µ b (see Figure 1). It is well known that this set of digital points denoted by L is a simple 8connected digital straight line (DSL) that can be defined by the diophantine equation [START_REF] Reveilles | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF][START_REF] Debled-Rennesson | A linear algorithm for segmentation of digital curves[END_REF]

0 ≤ ax -by + µ < max (|a|, |b|) = b.
The slope of L is a b , µ is the shift at origin. The remainder of L for a given digital point (x, y) is the value axby + µ.

Two particular straight lines are associated to this DSL L:

• the upper leaning line given by equation axby + µ = 0 and,

• the lower leaning line given by equation axby + µ = b -1.

The digital points lying on these lines are similarly called leaning points (see Figure 2). Upper leaning points have remainder value 0 while lower leaning points have remainder value b -1 (see [START_REF] Coeurjolly | Géométrie discrète et images numériques[END_REF] for more details). Since a, b, µ are integers and L is digitized with OBQ, L is identified as the upper leaning line of L.

3 Optimal covering of a straight line

Setting the problem

The objective is to cover the straight line L with closed balls centered on the points of the DSL L. Moreover, the union of these balls shall not contain any other digital point. These balls are defined according to the L ∞ metric such that the ball of radius ε is defined by B(P, ε) = {q|d ∞ (P, q) ≤ ε}. In our framework, we assume that the radius ε is the same for all the points of L (see Figure 3). In the general case, setting ε to 1 2 is not sufficient to cover the straight line L except for the special cases of horizontal, vertical and diagonal lines. We can notice the following elementary property: Property 1. The band delimitated by the 2 leaning lines (upper and lower) has a vertical thickness of b-1 b . The proof is straightforward from the equations of the leaning lines. Figure 3 illustrates the covering of a straight line of slope a b by balls of radius b-1 b . We remark that this radius is not sufficient since parts of the straight line remain uncovered.

To analyze this covering, we proceed by successive couples of adjacent pixels. It is easy to verify that the vertical distance between a pixel of remainder r of the DSL and the line L is equal to r/b, with r varying from 0 to b -1. The value b-1 b is obtained for a lower leaning point. Let us consider two successive pixels of L such that the first one is a lower leaning point. Then we notice that the union of the balls of radius b-1 b centered on these points leaves a part of L uncovered (see Figure 3). So the minimum ε for the line L to be covered is greater than b-1 b . In the following, we denote by P 0 , P 1 , . . . , P n the ordered sequence (increasing abscissae) of pixels of the DSL L. For each pixel P i we can define as proj(P i ) the vertical projection of P i on the straight line L. Consequently, the Euclidean straight line L may be partitionned as the union of subsegments [proj(P i ), proj(P i+1 )].

Property 2. If for every pair of successive pixels P i and P i+1 , B(P i , ε) ∪ B(P i+1 , ε) covers the straight segment [proj(P i ), proj(P i+1 )] then the straight line L is covered by the union of the balls B(P i , ε) .

As said before, ε must be greater than b-1 b , and it is easy to see that it must also be strictly lower than 1, otherwise new digital points are included in the union of balls.

U 2 L 2 L 1 b-1 b

Optimal covering

The following theorem defines the optimal value of ε as a function of the straight line parameters a and b.

Theorem 1. Let L a straight line of equation axby + µ = 0 and L its digitization with the OBQ scheme. In the context of the previous definitions, the union of balls B(P i , ε) centered on pixels of the DSL L with radius ε = max ( 1 2 , |a|+|b|-1 |a|+|b| ) covers the straight line L. This set doesn't contain any other digital pixels excepted those of the DSL.

Proof. We suppose that 0 ≤ a ≤ b. First of all, if the parameters of L are (0, 1, µ) (horizontal straight line) or (1, 1, µ) (diagonal straight line), the optimal value of ε is trivially equal to 1 2 . Otherwise, b is greater than or equal to 2. The distance between a point P i and its projection proj(P i ) (see above) is equal to r b if the remainder of P i is equal to r. In the case "P i and P i+1 4-connected", if r is the remainder of P i then r + a is the remainder of P i+1 . Consequently, the distance between P i+1 and its projection proj(P i+1 ) is equal to r+a b . We have the following inequalities:

• b-1 b > r b (P i belongs to L), • b-1 b ≥ r+a b (P i+1 belongs to L), • 2 b-1 b ≥ 1 (since b ≥ 2)
. With the first two inequalities, we deduce that proj(P i ) and proj(P ) covers the segment [proj(P i ), proj(P i+1 )] (see Figure 4). In the case "P i and P i+1 8-connected", according to the Figure 5(b), we introduce the point Q belonging to the straight line L situed at equal distance (L ∞ norm) from P i and P i+1 . The point Q is on L and shall belong to the two balls B(P i , ε) and B(P i+1 , ε). If r is the remainder of P i , we denote by ε(r) the minimum radius such that B(P i , ε(r)) ∪ B(P i+1 , ε(r)) covers the segment [proj(P i ), proj(P i+1 )]. We can write ε(r) as the sum r b + ρ(r), with ρ(r) ≥ 0 (see Figure 5). Using Thalès theorem, we have

ρ(r) a(b-r) b 2 = b-r b -ρ(r) b-r b
We obtain

ρ(r) = a(b -r) b(a + b)
Then ε(r) = ρ(r) + r b = a+r a+b . ε(r) is increasing and is maximum for r = b-1, corresponding to ε = a+b-1 a+b . This proof was developed in the case 0 ≤ a ≤ b. It can be easily extended in the general case to obtain ε = max( 1 2 , |a|+|b|-1 |a|+|b| ). Since ε < 1, P i is the only digital point in the ball B(P i , ε) for all i, which ends the proof.

Discrete Convexity

Definitions

An important literature has been developed for digital convexity [START_REF] Sklansky | Recognition of convex blobs[END_REF][START_REF] Kim | Digital and cellular convexity[END_REF][START_REF] Kim | Digital convexity, straightness, and convex polygons[END_REF][START_REF] Kim | Digital straight lines and convexity of digital regions[END_REF][START_REF] Chassery | Discrete convexity: Definition, parametrization, and compatibility with continuous convexity[END_REF]. A large number of these definitions are considered to be equivalent in case of simply connected sets [START_REF] Ronse | A bibliography on digital and computational convexity[END_REF][START_REF] Eckhardt | Digital lines and digital convexity[END_REF].The following definition is issued from the transcription of convexity definition in Euclidean space to discrete space. In Euclidean geometry, a region R is convex if and only if for every pair of points p, q belonging to R, the straight line segment [p, q] is included in R. The following definition of discrete convexity replaces inclusion by covering with balls [START_REF] Chassery | Discrete convexity: Definition, parametrization, and compatibility with continuous convexity[END_REF].

Definition 1 (ε-convexity). A connected component S is ε-convex, with ε be- longing to interval [ 1 2 , 1[ if:
• for every pair of pixels P and Q of S,

• for every real value α belonging to ]0, 1[, there exists a pixel R belonging to S such that the point (αP + (1α)Q) of the straight line supported by the two points P Q belongs to the balls B(R, ε), centered on R with radius ε.

Definition 2 (Discrete convexity)

. A connected component S is discrete convex if there exists a real ε ∈ [ 1 2 , 1[ such that S is ε-convex.

Algorithmic approach

The computation of the convex hull Conv(S) is done in two steps :

1. first, the xy convex shape issued from S is computed;

2. then the convex hull is computed from the xy convex shape.

The xy convex shape of a connected component is defined as the convex shape along horizontal and vertical directions: for any points P and Q of the xy convex shape, if P and Q are on the same line or column of the grid, then all the points between P and Q (on this line or column) also belong to the xy convex shape. It is evident that the xy convex shape is included in the convex hull.

Algorithm 1 describes how to compute the x -y convex shape. Suppose that the shape S is included in a binary image of size m × n. The initialization step consists in: 1. compute the indices i f and i l of the upper and lower lines containing points of S respectively; 2. for each line i between i f and i l , compute min(i) and max(i) as the minimum and maximmum indices of points of S on line i.

This initilization step is done in O(nm). Afterwards, the three loops of Algorithm 1 are performed to compute the maximal and minimal abscissas for each line of the xy convex shape of S.

Algorithm 1: x-y convex shape computation of a set S of discrete points

for i from i f to i l do 1 min ↓ (i) = min(min(i), min ↓ (i -1)) ; // min ↓ (i) is initially set to m max ↓ (i) = max(max(i), max ↓ (i -1)) ; // max ↓ (i) is initially set to 0 end for i from i l to i f do 2 min ↑ (i) = min(min(i), min ↑ (i + 1)) ; // min ↑ (i) is initially set to m max ↑ (i) = max(max(i), max ↑ (i + 1)) ; // max ↑ (i) is initially set to 0 end for i from i f to i l do 3 min(i) = max(min ↓ (i), min ↑ (i)) max(i) = min(max ↓ (i), max ↑ (i)) j f = min(j f , min(i)) ;
// j f is initially set to n j l = max(j l , max(i)) ; // j l is initially set to 0 end At the end of Algorithm 1, the variables min(i) and max(i) contain the indices of the minimum and maximum points of the xy convex shape for the line i. This algorithm also computes the coordinates i f , i l , j f , j l of the bounding rectangle of S.

Finally, the convex hull of S is computed as follows. To do so, we define the extremal point of the bounding rectangle as the points of S belonging to this rectangle and such that at least one of its two neighbours on the rectangle does not belong to S. These eight points are depicted on Figure 6(d).

Finally, the convex hull is computed from these extremal points using a technique similar to the algorithm of Sklansky [START_REF] Sklansky | Measuring concavity on a rectangular mosaic[END_REF][8](Chap. 13) on the polygon defined by the variables min(.) and max(.). Note that this polygon is simple and completely visible from the outside, such that the algorithm works in this case. This algorithm works in O(max(m, n)) time, which leads to a global complexity of O(nm) to compute the convex hull of S. The result on the example of Figure 6 is given in Figure 7.

Optimal convering and discrete convexity

Let us denote by {P i , i = 0..n} the ordered set of vertices of the convex hull Conv(S). To each edge [P i , P i + 1] we associate the slope parameters (a i , b i ) such that gcd(a i , b i ) = 1. From Property 1, this edge is covered by balls centered on the pixels of its OBQ digitization with radius max( 1 2 , |ai|+|bi|-1 |ai|+|bi| ).

It is easy to prove that the connected set Conv(S) of discrete points included in Conv(S) is ε-convex with ε = max i∈0...n max( Moreover such definition is fully compatible with the continuous one as it is proved by the theorem [1]: Theorem 2. Let S be a connected component in space R 2 . If, for every sampling step, the discrete connected component S attached to S is discrete convex, then S is convex in space R 2 .

Using S and Conv(S) a lot of features can be extracted to estimate a measure of convexity of S (number of missing pixels, distribution of concavities) [START_REF] Roussillon | Algorithmes d'extraction de modèles géométriques discrets pour la représentation robuste des formes[END_REF]. These features are used to measure the degree of convexity of a shape in presence of many concavities.

Conclusion

The transcription of continuous concept like convexity into digital space needs specific attention for its definition as well as for its properties. In case of convexity we use the definition of digital straight line and we replace the notion of inclusion by the notion of covering. The main result obtained in this paper is issued from parametric representation of digital straight line to characterize the optimal radius for covering balls centered on digital straight line pixels.

D

Figure 1 :

 1 Figure 1: A straight line L of equation 3x -8y + µ = 0 and its OBQ digitization analytically given by the digital straight line of equation 0 ≤ 3x -8y + µ < 8.

Figure 2 :

 2 Figure 2: Digital straight line of equation 0 ≤ 3x -8y + µ < 8 with upper (lower) leaning points U i (L i ).

Figure 3 :

 3 Figure 3: Balls of radius ε = b-1 b centered on the DSL points do not cover L: uncovered regions are indicated by arrows.

  i+1 ) are covered by B(P i , b-1 b ) ∪ B(P i+1 , b-1 b ). The third inequality ensures that B(P i , b-1 b ) and B(P i+1 , b-1 b ) overlap. We can conclude that the union of the two balls B(P i , b-1 b ) ∪ B(P i+1 , b-1 b

Figure 4 :

 4 Figure 4: When P i and P i+1 are 4-connected, B(P i , b-1 b ) ∪ B(P i+1 , b-1 b ) covers the segment [proj(P i ), proj(P i+1 )].

Figure 5 :

 5 Figure 5: Illustration of the proof when P i and P i+1 are 8-connected. r is the remainder of P i . (b) is a close up of the region circled in (a).

Figure 6 (

 6 Figure 6(b) illustrates the first two loops of the algorithm (lines 1 and 2).Figure 6(c) illustrates the third loop (line 3).At the end of Algorithm 1, the variables min(i) and max(i) contain the indices of the minimum and maximum points of the xy convex shape for the line i. This algorithm also computes the coordinates i f , i l , j f , j l of the bounding rectangle of S.Finally, the convex hull of S is computed as follows. To do so, we define the extremal point of the bounding rectangle as the points of S belonging to this rectangle and such that at least one of its two neighbours on the rectangle does not belong to S. These eight points are depicted on Figure6(d).

Figure 6 (

 6 Figure 6(b) illustrates the first two loops of the algorithm (lines 1 and 2).Figure 6(c) illustrates the third loop (line 3).At the end of Algorithm 1, the variables min(i) and max(i) contain the indices of the minimum and maximum points of the xy convex shape for the line i. This algorithm also computes the coordinates i f , i l , j f , j l of the bounding rectangle of S.Finally, the convex hull of S is computed as follows. To do so, we define the extremal point of the bounding rectangle as the points of S belonging to this rectangle and such that at least one of its two neighbours on the rectangle does not belong to S. These eight points are depicted on Figure6(d).

Figure 6 :

 6 Figure 6: (a) A discrete shape S. (b) Illustration of the variables min ↓ , min ↑ ,max ↓ , max ↑ used in Algorithm 1. (c) Illustration of the variables min(.) and max(.) at the end of Algorithm 1. (d) Extremal points of S.

Figure 7 :

 7 Figure 7: Convex hull of the discrete shape S: the point marked by a square is added during the application of Sklansky's algorithm. The discrete points belonging to the polygonal line are the vertices of Conv(S). Conv(S) is defined as the discrete points inside Conv(S).

  1 2 , |ai|+|bi|-1 |ai|+|bi| ). S is discrete convex if and only if it is equal to Conv(S).