Laurent Desbat
email: laurent.desbat@imag.fr

Simon Rit
email: simon.rit@gmail.com

R Clackdoyle

C Mennessier

Emmanuel Promayon

S Ntalampeki

Algebraic and analytic reconstruction methods for dynamic tomography

Algebraic and analytic reconstruction methods for dynamic tomography.

I. INTRODUCTION

Dynamic tomography is a very active area [START_REF] Grass | Helical cardiac cone beam reconstruction using retrospective ECG gating[END_REF], [START_REF] Ford | Respiratorycorrelated spiral CT: a method of measuring respiratory-induced anatomic motion for radiation treatment planning[END_REF]. Movements or deformations of the patient must be taken into account in nuclear medicine, such as SPECT or PET because of long measurement time, but also in CT for fast moving organs such as the heart. Generally, patient movements or deformations occur in 3D. In this work we consider dynamic 3D Cone Beam tomography. We present both algebraic and analytic aproaches to compensate for the patient deformation during the reconstruction.

A. 3D cone beam notations

Let f : R 3 → R be the 3D attenuation distribution to be reconstructed from projections (x-ray, SPECT or PET). In this work, we consider mainly the 3D cone beam transform [START_REF] Grass | Helical cardiac cone beam reconstruction using retrospective ECG gating[END_REF] where ζ ∈ S 2 is a unit vector in R 3 (S 2 is the unit sphere in R 3 ), a(t) ∈ R 3 is the x-ray source position at time t ∈ T ⊂ R, t is then also the source trajectory parameter in R 3 , see Fig. 1. D t f ( ζ), at fixed t, is assumed to be acquired in a negligible time, ∀ ζ ∈ S 2 . The function f is assumed to have a compact support and a(t) has a strictly positive distance to the support of f . This transform appears in 3D x-ray tomography (reconstruction from 2D xray projections or multiline CT) with applications in cardiac CT or radiotherapy.

g D (t, ζ) def = Df (t, ζ) def = D t f ( ζ) def = +∞ 0 f a(t) + l ζ dl,
3D cone beam tomography problem is the reconstruction of f from g D . These last years, new developments have been proposed to solve analytically, exactly, and efficiently this problem, in particular for the helical source trajectory, but also for more general trajectory (see for exemple [START_REF] Grangeat | Mathematical framework of cone beam 3D reconstruction via the first derivative of the Radon transform[END_REF], [START_REF] Defrise | The solution to the 3D image reconstruction problem from 2D parallel projections[END_REF], [START_REF] Katsevich | Analysis of an exact inversion algorithm for spiral cone beam CT[END_REF], [START_REF]Theoretically exact filtered back-projection type inversion algorithm for spiral CT[END_REF], [START_REF] Pack | Cone-beam reconstruction using 1d filtering along the projection of m-lines[END_REF], [START_REF] Zou | Image reconstruction on PI-lines by use of filtered backprojection in helical cone-beam CT[END_REF]).

rolf.clackdoyle@univ-st-etienne.fr, mennessier@cpe.fr

B. Deformation model

In dynamic tomography, the attenuation f to be reconstructed is also a function of t. In our approach, we consider a time dependent deformation of the space as in [START_REF] Crawford | Respiratory compensation in projection imaging using a magnification and displacement model[END_REF], [START_REF] Roux | Exact reconstruction in 2D dynamic CT: compensation of time-dependent affine deformations[END_REF]. We introduce a time dependent deformation model Γ t and we assume that, for t ∈ T , Γ t are known bijective appropriately smooth functions on R 3 whose inverse are smooth too. We assume that the attenuation at point x ∈ R 3 at time t can be written in the form f Γ t ( x) = f Γ t ( x) , where f is the attenuation reference point (for example, at t = 0). So Γ t simply transforms the position x at time t to its position Γ t ( x) at the reference time (t = 0).

The problem under consideration is the reconstruction of f , or equivalently of f Γ t , from cone-beam measurements Df Γt , see [START_REF] Grass | Helical cardiac cone beam reconstruction using retrospective ECG gating[END_REF], for known deformation functions Γ t , t ∈ T . Our objective is to apply efficient accurate deformation corrections in the reconstruction algorithms.

C. Dynamic reconstruction approaches

General patient deformations can be nowadays compensated only by algebraic reconstruction approaches. Iterative reconstruction requires an adequate projector able to perform the attenuation sum (1) based on a discrete representation of f at the reference time. This is done in the static case by computing the intersection of the straight acquisition lines with basis functions, e.g. voxels. In the dynamic case, the projector must compensate for the motion. If Γ t preserves the rectitude of acquisition lines, i.e. if it transforms a set of convergent aquisition lines at time t into a set of convergent lines at the reference time, this can be done as in the static case along the virtual transformed lines. In particular, if the patient motion is globally rigid, this amounts to using virtual source and detector positions and applying the reconstruction algorithm as in the static case [START_REF] Fulton | Use of 3D reconstruction to correct for patient motion in SPECT[END_REF], [START_REF] Hutton | A hybrid 3-D reconstruction/registration algorithm for correction of head motion in emission tomography[END_REF], [START_REF] Feng | Use of three-dimensional Gaussian interpolation in the projector/backprojector pair of iterative reconstruction for compensation of known rigid-body motion in SPECT[END_REF] or to deform data prior to reconstruction [START_REF] Livieratos | Rigid-body transformation of list-mode projection data for respiratory motion correction in cardiac PET[END_REF], [START_REF] Lamare | Respiratory motion correction for PET oncology applications using affine transformation of list mode data[END_REF]. These approaches can be generalized to deformations transforming the set of acquisition lines of each cone beam projection into other sets of concurent lines [START_REF] Crawford | Respiratory compensation in projection imaging using a magnification and displacement model[END_REF], [START_REF] Roux | Exact reconstruction in 2D dynamic CT: compensation of time-dependent affine deformations[END_REF], [START_REF] Desbat | Compensation of some time dependent deformations in tomography[END_REF]. Thus, analytic approaches essentially allow for the compensation of deformations in subclasses of those presented in section I-B. With this approach, deformations leaving globally invariant the cone beam geometry acquisition can be compensated. When this is not the case, iterative reconstruction is still feasible but the intersection of the straight acquisition line with f Γt ( x) should take into account the translation and the deformation of the basis functions. Moreover, in 3D interventional image reconstruction from x-ray projections, the source trajectory is very often a circle [START_REF] Desbat | Minimally Invasive Interventional Imaging for Computer Assisted Orthopedic Surgery[END_REF], [START_REF] Wiesent | Enhanced 3-d-reconstruction algorithm for c-arm systems suitable for interventional procedures[END_REF] for which exact analytic inversion does not apply whereas algebraic approaches are well defined. When a high contrast object like coronary arteries is being reconstructed, the deformation can be ignored [START_REF] Blondel | 3D tomographic reconstruction of coronary arteries using a precomputed 4D motion field[END_REF], [START_REF] Blondel | Reconstruction of coronary arteries from a single rotational X-ray projection sequence[END_REF]. Otherwise, one has to find a way to deform the basis functions. For example, spherical basis functions can be deformed in ellipoids as a first approximation [START_REF] Reyes | Respiratory motion correction in emission tomography image reconstruction[END_REF], [START_REF] Reyes | Respiratory motion correction in emission tomography imaging[END_REF].

In the next section, we present a dynamic tomography algebraic method framework. In section III, we present an analytic approach for dynamic tomography. Then, we present some numerical experiments in section IV.

II. ALGEBRAIC DYNAMIC RECONSTRUCTION

In algebraic approch, we assume that the unknown function f can decomposed into a finite linear combination of given independent functions, such as voxel indicator functions. Let (e j ) j∈J be a set of basis functions e j : R 3 → R, J being a finite set of index such as J = {1, . . . , n J }, n J ∈ N denoting the number of elements #(J) of the set J (multi index sets are also very often used in multi dimension space). For example, e j can be the indicator function of the voxel number j in the reconstruction region. We assume that

f ( x) = j∈J f j e j ( x), (2)
where f j ∈ R, j ∈ J. We also assume that the acquisition can be modeled by

d i = Ω h i ( x)f ( x)d x, i ∈ I. ( 3
)
where (d i ) i∈I is the real vector of acquired data. In 2D tomography, h i ( x) could be the dirac on a line δ( x

• θ i -s i )
where θ i is the direction of the projection and s i ∈ R is the signed distance of the line to the center. In 3D h i ( x) could be the dirac on the x-ray line but it could also be the indicator of the conical region joining a point source a i and a detector or some smooth response function obtained by calibration of the X-ray, PET or SPECT system. In nuclear imaging, h i can also model more physics such as attenuation, etc. The number of data n I is finite.

In static tomography, the algebraic aproach yields a linear system to be solved

d i = Ω h i ( x)f ( x)d x = Ω h i ( x) j∈J f j e j ( x)d x = j∈J Ω h i ( x)e j ( x)d x f j (4)
We have to solve the linear system d = Af where d = (d i ) i=1,...,n I is the known vector of data, f = (f j ) j=1,...,n J is the unknown vector of coefficients of f to be identified and the matrix entry A i,j is A i,j = Ω h i ( x)e j ( x)d x. The matrix is generally sparse because both h i and e j functions have a limited support in the domain Ω ⊂ R 3 . Thus iterative methods are used to solve the linear system.

In dynamic tomography

d t,i = Ω h i ( x)f Γ t ( x) d x = j∈J Ω h i ( x)e j Γ t ( x) d x f j ( 5
)
Let us now assume that the functions e j Γ t ( x) can be decomposed (or approximated) into a finite linear combination of given independent functions (b k ) k∈K , K being a finite index (or multi-index) set, more precisely

e j Γ t ( x) = k∈K Γ t k,j b k ( x). (6)
then [START_REF]Theoretically exact filtered back-projection type inversion algorithm for spiral CT[END_REF] introduced in (5) yields

d t,i = j∈J Ω h i ( x) k∈K Γ t k,j b k ( x)d x f j = j∈J k∈K B i,k Γ t k,j f j ( 7
)
where B is the algebraic matrix for the basis function

(b k ) k∈K B i,k = Ω h i ( x)b k ( x)d x. (8)
In 3D CB tomography, i is usually a multi-index: at least one index, say i 1 , is related to the source position on its trajectory, and an other index (or multi-index) i 2 is related to the detector pixel position in space at the source position i 1 or equivalently the direction of the ray from the source position at t(i 1 ) (usually, in 3D CB, i 2 is a multi-index of two values because the direction space S 2 is a two dimensional set). Thus t is a function of i 1 (the time depends on the source position), thus d t,i is d t(i1),i2 , where i

= (i 1 , i 2 ), i 1 ∈ I 1 , i 2 ∈ I 2 , I = I 1 × I 2 and d t(i1),i2 = Df Γ t(i 1 ) t(i 1 ), ζ(i 2 ) = j∈J k∈K B (i 1 ,i 2 ),k Γ t(i 1 ) k,j f j (9)
Thus,

d i 1 = B i 1 Γ t(i 1 ) f , i 1 = 1, . . . , n I 1 ( 10
)
where the n I 1 matrices B i 1 are n I 2 × n K (n I 1 sets of n I 2 lines of a classical ART matrix B corresponding to each i 1 source position, i 1 = 1, . . . , n I 1 ) and the matrices Γ t(i1) are n K × n J matrices coding the deformations by Γ t(i 1 ) of the basis functions e j into the basis fuctions b k , for each considered time t(i 1 ), i 1 = 1, . . . , n I 1 . With this formulation, an algebraic matrix for dynamic tomography differs from the classical algebraic method only by the introduction of the matrices Γ t(i 1 ) . Just as for classical tomography, iterative methods are used for solving the linear system [START_REF] Roux | Exact reconstruction in 2D dynamic CT: compensation of time-dependent affine deformations[END_REF] and the matrices B and Γ t(i1) , i 1 = 1, . . . , n I 1 are generally not stored but computed during the iterations.

III. ANALYTIC DYNAMIC RECONSTRUCTION

We recall briefly in this section a class of 3D deformations which can be analytically compensated in 3D Cone Beam reconstruction, see [START_REF] Desbat | Compensation de déformations en tomographie dynamique 3D conique[END_REF] and [START_REF] Desbat | Compensation of some time dependent deformations in tomography[END_REF], [START_REF] Roux | Exact reconstruction in 2D dynamic CT: compensation of time-dependent affine deformations[END_REF] for 2D fan-beam dynamic tomography. In order to stay in the framework of 3D CB reconstruction, we consider the deformations Γ t which transform the set of convergent acquisition lines at time t into an other set of convergent lines at the reference time, i.e. we consider deformations that globally preserve the CB acquisition geometry. The deformation Γ t : R 3 → R 3 preserves the CB aquisition geometry globally if it transforms the source point a(t) ∈ R 3 at time t into a virtual source point Γ t ( a(t)) at the reference time (e.g. t = 0 as we assumed) and if any half line from a(t) is transformed into a virtual half line from the virtual source Γ t ( a(t)), see Fig. 2. Let us use the spherical coordinates (l, ζ) ∈ R + × S 2 of x -a(t), i.e., x = a(t) + l ζ. The half line from a(t) in the direction ζ is denoted by a(t

) + R + ζ. Assume x belongs to a(t) + R + ζ. Γ t leaves the CB geometry globally invariant if ∀t ∈ T, ∀ ζ ∈ S 2 , Γ t a(t) + R + ζ = Γ t ( a(t)) + R + Γ S 2 ,t ( ζ), ( 11
)
where Γ S 2 ,t S 2 -→ S 2 is a diffeomorphism (bi-regular bijection) on the unit sphere which associates to a direction ζ at t a direction Γ S 2 ,t ( ζ) at the reference time. More precisely, the following deformation leaves the 3D CB geometry globally invariant:

Γ t ( x) = Γ t a(t) + l ζ = Γ t ( a(t)) + Γ t, ζ (l) Γ S 2 ,t ( ζ), (12)
where Γ t, ζ is a bi-regular bijective function on R + such that Γ t, ζ (0) = 0. In the following Γ t, ζ is linear in order to stay in the framework of the 3D CB transform (more complex bijections would yield a non constant Jacobian in [START_REF] Desbat | Compensation of some time dependent deformations in tomography[END_REF] which would lead to generalized CB transforms for which we do not have inversion formulas). Thus let Γ t, ζ (l) = c t, ζ l with c t, ζ > 0 being both a function of t but also of ζ. We then can write [START_REF] Hutton | A hybrid 3-D reconstruction/registration algorithm for correction of head motion in emission tomography[END_REF] as:

Γ t ( x) = Γ t ( a(t)) + c t, ζ l Γ S 2 ,t ( ζ) (13)
Now, let v(t) def = Γ t ( a(t)) -a(t). We can decompose Γ t of (13) into:

Γ t = T v(t) • ∆ t , ( 14
) x 2 x 3 ζ x x 1 l x 2 x 1 x 3 Γ t Γ t ( x) Γ S 2 ,t ( ζ) Γ t ( a(t)) v(t) a(t) a(t) c t, ζ l
Fig. 2. Deformation Γ t preserving the 3D CB geometry globally where T v(t) is the translation of vector v(t)

( T v(t) ( x) = v(t) + x) and ∆ t ( x) = a(t) + c t, ζ l Γ S 2 ,t ( ζ). ( 15
)
We then remark that the deformation ∆ t can be analytically compensated directly within each 3D projection at fixed t:

D t f ∆t ( ζ) = +∞ 0 f a(t) + c t, ζ l Γ S 2 ,t ( ζ) dl = 1 c t, ζ D t f ( Γ S 2 ,t ( ζ)). ( 16
)
Thus, from [START_REF] Desbat | Compensation of some time dependent deformations in tomography[END_REF] we have

D t f ( ζ) = c t, Γ -1 S 2 ,t ( ζ) D t f ∆ t ( Γ -1 S 2 ,t ( ζ)). ( 17
)
Combining now ( 14) and ( 17) we have

D t f T v(t) ( ζ) = c t, Γ -1 S 2 ,t ( ζ) D t f T v(t) • ∆t ( Γ -1 S 2 ,t ( ζ)). ( 18
)
The deformation ∆ t can thus be very simply compensated within the projection

D t f T v(t)
D t f T v(t) ( ζ) = +∞ 0 f a(t) + v(t) + l ζ dl = +∞ 0 f Γ t ( a(t)) + l ζ dl. ( 19
)
Thus, the Tuy-Grangeat stable reconstruction conditions [START_REF] Tuy | An inversion formula for cone-beam reconstruction[END_REF],

[25] at a point x will be read on the virtual trajectory Γ t ( a(t)), just as the possibility to use modern 3D CB reconstruction algorithms [START_REF] Katsevich | A general scheme for contructing inversion algorithms for cone beam CT[END_REF], [START_REF] Pack | Cone-beam reconstruction using 1d filtering along the projection of m-lines[END_REF], [START_REF] Sidky | A minimum data FBP-type algorithm for image reconstruction in cone-beam CT[END_REF].

IV. NUMERICAL EXPERIMENTS

A realistic digital phantom of the thorax was used to evaluate the behavior of different reconstruction methods in the presence of respiratory motion. Lung radiotherapy optimization is a typical application for this study. The phantom was obtained from a 4D CT image acquired on a scanner synchronised with a respiratory signal [START_REF] Vandemeulebroucke | Point-validated pixel-based breathing thorax model[END_REF]. One of its ten 3D CT images was chosen to represent f . Nine dense motion vector fields were computed by deformable registration between each 3D CT image and the reference. From them, continuous trajectories for each voxel of f were obtained by supposing that the motion is linear between two consecutive respiratory instants. A set of cone-beam projections g D were computed using a projector taking into account this motion model and the geometry of an existing cone-beam CT scanner. Fig. 3 allows the comparison of the reference CT image with different reconstructions.

Motion artifacts, such as blur and streaks, are clearly visible, both with Feldkamp's analytic method [START_REF] Feldkamp | Practical cone-beam algorithm[END_REF] and with Simultaneous Algebraic Reconstruction Technique (SART) [30], which suppose that the patient remained static during the acquisition. The motion does not globally preserve the cone-beam geometry as described in part III, that is why only a heuristic motion compensation during the backprojection step of Feldkamp's algorithm was feasible [31]. The blur is reduced but streaks are still visible, particularly on the axial slice, because the compensation is not exact. Finally, motion compensated SART with the projector used for the computation of the cone-beam projections fully eliminate the motion artifacts.

V. CONCLUSIONS AND FUTURE WORKS

We have presented a framework for both analytic and algebraic method in dynamic tomography with known time dependent deformations of the space. We have provided numerical experiments showing that algebraic methods can compensate the respiratory motion in 4D CT from acquisition during lung radiotherapy applications. In future works, we will compare analytic and algebraic approaches on this application.

Fig. 1 .

 1 Fig. 1. Cone Beam geometry parameters

  • ∆ t , i.e. within D t f Γ t , in order to compute D t f T v(t) . Afterwards the translation T v(t) can be simply compensated: instead of reconstructing f from the 3D CB acquisition on the real trajectory a(t), t ∈ T ⊂ R, we have to perform the reconstruction from the virtual trajectory Γ

t ( a(t)) = a(t) + v(t), t ∈ T ⊂ R. Indeed, from the definition

[START_REF] Grass | Helical cardiac cone beam reconstruction using retrospective ECG gating[END_REF]

This work was supported by the ANR grant NT05-1_45428, ToRIID


Alexandre Carra

Jean-Louis Martiel

Emmanuel Promayon

Virtual adhesive: a way to handle sticky collisions in surgical and biological simulators

Keywords: (according to ACM CCS): I.3.5 [Computer Graphics]: Physically based modeling I.3.7 [Three-Dimensional Graphics and Realism]: Animation Physically-based model, sticky contact modeling, collision response, collision detection, deformable objects

A variety of methods have been proposed to efficiently process collisions between deformable objects. The method presented in this paper allows to model sticky states between deformable objects with triangulated surfaces. In contrast to an often used approach that consists of generating forces which eventually separate colliding objects, our method is based on the creation of an adhesive virtual object (virtual adhesive). This virtual adhesive is composed by "clones" of all particles locally in collision. Particle clones are used to gather forces on the virtual adhesive, which behaves like a rigid body. The resulting displacement of the virtual adhesive is used in turn to constrain the particles displacement: particles stick to their clone. As a result, no further interpenetration is possible and a sticky state is obtained in the considered zone. Moreover, the method correctly resolves contacts without introducing additional energy in the system. Results are presented through several simulations.

Introduction

Collision detection and response between deformable bodies have been studied quite extensively [TKH * 05]. They are crucial elements of surgical simulators, especially for modeling the interactions between organs and surgical tools, such as lancets and needles. Biological and cellular simulators, also called virtual cell modeling and simulation, are another type of emerging tools in which the collision detection and response is required to model cell deformations and migrations [START_REF] Promayon | Physically-based 3d simulations of cell deformations and migrations[END_REF]. Although the objectives of surgical simulators are very well understood today, their biological counterpart is of a different nature.

Cell movements and deformations are essential for many important biological and medical processes, including cell division and migration, the formation of viable embryos or the growth of cancerous tumours. A long-standing challenge in modeling complex biological systems is to combine mechanisms of different physical origins (e.g. cell elasticity, cytoskeleton dynamics, chemical reactions) at different times (from the second to the minute) and length scales (from micrometers to tens of millimetres). If we adopt a reductionist approach, the force generation in cells results from the dynamic assembly of individual molecules or filaments (actin, tubulin). These physico-chemical processes are now well understood [START_REF] Pollard | Cellular motility driven by assembly and disassembly of actin filaments[END_REF]. For example, actin filaments produce piconewton forces, as shown in experiments [START_REF] Kovar | Insertional assembly of actin filament barbed ends in association with formins produces piconewton forces[END_REF] or in physical modeling [BMB * 07]. Although these models are necessary to decipher the resulting behavior of complex macromolecular assemblies, they are by far too complex to be useful to account for the dynamics of a cell or a cell population. Therefore, the next step in simplification of approach would use phenomenological models, based on a global description of phenomena, without relying on precise and detailed mechanisms in cells. For example, interactions between a cell rolling on a surface by creating/destroying c The Eurographics Association 2007. bonds between the membrane and proteins on the surface can be handled by kinetic models [START_REF] Dembo | The reaction-limited kinetics of membrane-to-surface adhesion and detachment. Proceedings of the Royal Society of London[END_REF]. These models are based on a phenomenological description of cells rolling over rigid surfaces or over each other and do no necessitate complex mathematical/computer solutions. Depending on the question addressed and on the required degree of accuracy of the representation of experimental data, we can use different approaches, from the simulation of partial differential equations (e.g. to get the elastic response of cells or tissues; to represent the different levels of molecules in cells) to the discrete iterations of an algorithm that finds the problem solution by minimizing a global energy function. However, we can go a step further in the simplification in the mathematical/computer-based modeling. The question "How do cells move on solid substratum or interact with other cells?" can also be analyzed via very simple computer models. In addition, these approaches give valuable insights into the global cellular mechanism that permits cells to deform, slide roll or stall, when they interact with other cells or rigid objects.

In both surgical and biological simulators, contact must be handled efficiently. However a category of contact that are very rarely considered in simulators, due to the complexity of the phenomenon and because there exists no simple theoretical model in physics today is sticky contacts. Deformable organs, tissues or cells are generally surrounded by a membrane or capsule that have adhesive properties. These adhesive properties considerably modify the contact interactions. In such simulators, contact generally occurs for quite a long time and forces between the colliding bodies have to be exchanged. Furthermore, to model adhesive properties a kind of cohesion between the colliding bodies also have to be expressed, which has to be maintained as long as the relative deformation and movement of the bodies and their intrinsic membrane adhesive properties enable it. This paper presents an alternative solution to modeling sticking deformable objects. Based on a phenomenological approach, it consists of creating a virtual adhesive body in the zone where the bodies are colliding, in order to model the sticky contact. In our approach, a sticking zone is a zone where deformable objects are temporarily glued together i.e. stay in contact whitout intersecting while exchanging forces. Our goal is not to be physically or mechanically correct but instead to model sliding and sticking surfaces so that visual realism is maintained as well as good force transfer within a reasonable computation time.

Related work

The handling of contact for rigid and deformable bodies has been widely studied in computer graphics. In this work, we focus on collision between deformable objects. For an extended discussion, the readers can refer to the detailed paper on collision detection of deformable objects [TKH * 05].

Collision response schemes for deformable objects can be classified principally into two different categories: the penalty-based method and the constraints method. The penalty-based method is probably the most widely used because of its simplicity to implement. In this approach, a force is exerted on a colliding point whose magnitude is proportional to its penetration depth. Terzopoulos et al. has been one of the first to use it in [START_REF] Terzopoulos | Elastically deformable models[END_REF]. Later, Cani uses isopotential implicit surfaces generated by skeletons to model deformable solids [START_REF] Promayon | An implicit formulation for precise contact modeling between flexible solids[END_REF]. With this implicit formalism, exact contact surfaces are generated between objects. In [START_REF] Spillmann | Contact surface computation for coarsely sampled deformable objects[END_REF], Spillmann and Teschner also focus on the contact surface. They compute an explicit representation of the contact surface between flexible objects and triangulated surfaces. This contact surface exclusively depends on the internal forces; no user-defined parameters are needed. After collision response, the colliding vertices are displaced onto the contact surfaces and a collision-free configuration is obtained. With this method, the authors address discontinuity problems due to discrete time steps and coarse surface representations. Keiser et al. [KMH * 04] propose a combination of contact surface and penalty-based methods. They compute a virtual contact surface for colliding deformable objects and determine the position of the colliding points on the contact surface by applying penalty forces. In [HTK * 04], Heidelberger et al. present a method to calculate consistent penetration depth. They consider a set of close surface features to reduce discontinuities of the estimated penetration depth directions for small displacements of the penetrating vertices. The penetration depth and direction are only computed for the colliding points that are close to the penetration surface, and then propagated to the colliding points with larger penetrations. In [START_REF] Dequidt | Time-critical animation of deformable solids[END_REF], the authors present a multi-resolution model, based on a deformable octree structure that can handle large displacements. In this case, collision responses are computed using penalty forces. Recently, Barbic et al. [BJ07] have presented a model that can support haptic rendering of distributed contact between a rigid object and a reduced FEM deformable object (both with non-trivial geometry) by using a rigid distance field and a deformable pointshell. Penalty-based approach is employed to handle contacts: the penalty forces are determined by querying the points of the pointshell object against the signed distance field of the other object. The approach for time-critical evaluation of contact forces is based on several optimizations such as temporal coherence and graceful degradation.

Constraints can also be used to handle collision. The key idea is to define a system of constraints to prevent objects from penetration. This system is formulated as a Linear Complementary Problem (LCP) whose solution gives the contact forces. The contact forces are then computed by solving the system. Duriez et al. [START_REF] Duriez | Signorini's contact model for deformable objects in haptic simulations[END_REF][START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF] use Signorini's model to deal with contact between deformable objects. Then, after solving a LCP, the forces are integrated to find the constrained motion. Based on an explicit finite elc The Eurographics Association 2007. ement method, the approach developed in [DDCB01] locally adapts the sampling at which the simulation is computed to the current local deformation. When a collision is detected, the approach uses non-penetration constraints to move the collided surface points. Moreover, the active mesh nodes that are linked to a concerned surface point are also displaced. In [START_REF] Pai | Quasirigid objects in contact[END_REF], Pauly et al. combine the benefits of rigid body models for dynamic simulation and the advantages of deformable models for solving contacts and producing visible deformations. They set up a LCP that define the contact surface between two objects and solves the tractions that act on this surface. Then they compute the corresponding displacements and the total wrench that acts on the quasi-rigid bodies. In addition, volume preservation is ensured for any pressure distribution acting on the surface. Recently, in [START_REF] Spillmann | Non-iterative computation of contact forces for deformable objects[END_REF], the authors combine the constraints method with the penalty approach. They compute the positions of the colliding points and the contact surface by considering the stress on each point. A contact force is then calculated that accelerates the point immediately to its collision-free position. Here, the system of equations is decoupled and can be solved analytically. Moreover, the approach does not require the definition of response constants.

While all these approaches efficiently handle collision, they do not take into account sticky contact between deformable objects. One classic solution to this problem are approaches based on the coulomb friction law. Mirtich uses such a model to handle sticking and sliding contact between rigid bodies [START_REF] Mirtich | Impulse-Based Dynamic Simulation[END_REF][START_REF] Mirtich | Hybrid simulation: combining constraints and impulses[END_REF] but its solution is computationally expensive. In [START_REF] Clavet | Particlebased viscoelastic fluid simulation[END_REF], Clavet et al. present an impulsebased method to make particles stick. They supplement the impulse formulation with an attraction term that accounts for particles that are close to an object. Another solution consists of adding elastic links, e.g. springs between colliding nodes, and removing them when the surfaces move sufficiently appart, as was done in [START_REF] Promayon | Physically-based 3d simulations of cell deformations and migrations[END_REF]. However these three approaches introduce additional energy into the system.

Description of the model

This section first presents a general overview of the model used in this paper. The description of a rigid object is then detailed as it will be needed to model the virtual adhesive bodies. Finally, we explain the principle that allows to link a rigid and an elastic object together. Due to this constraint, an appropriate response is found when such two bodies are in contact.

General presentation

Our discrete model can handle objects which can be either elastic, rigid, or muscular [START_REF] Marchal | Simulating prostate surgical procedures with a discrete soft tissue model[END_REF][START_REF] Promayon | Physically-based 3d simulations of cell deformations and migrations[END_REF]. Each object is composed of a set of particles, which are defined by a position, a mass and a list of its neighbouring particles.

Neighbouring particles are initially defined by using a topological criterion that is invariant during the simulation. Particles are spread to volumetrically describe the object. The external surface of the objects are defined by a triangular mesh based on the surfacic particles which are the most external ones. Like in [BJ07], all the objects used in this paper are only described by their surfacic particles and mesh. Moreover the algorithms presented here work as well when particles are distributed in a volume. In the model, forces are applied on each particle. Two kinds of forces are considered here: force fields (e.g. gravitation force), and internal forces (e.g. elasticity).

The model also takes into account constraints associated to more complex properties such as the incompressibility of an object or the attachment between different bodies. Two kinds of constraints are implemented: first, local constraints can be applied to isolated particles. For instance, the position of a specific particle. Second, global constraints are used to implement incompressibility. For more details on the modeling of the elasticity and the method used to preserve the volume the readers can refer to [START_REF] Marchal | Simulating prostate surgical procedures with a discrete soft tissue model[END_REF].

The dynamics of the object is classically determined by the integration of the equations of motion. At each time step, boundary conditions are applied on the particles (e.g. null displacements). Next, the forces applied on each particle are computed and summed and the internal forces are redistributed to satisfy the action-reaction principle. The position of each particle is then determined by a standard explicit integration scheme. Finally, the local and global constraints are taken into account by adjusting these positions. The algorithm is detailed below for one iteration (Alg. 1).

Rigid object modeling

In this paragraph, the specific properties of a non deformable body are detailed. By definition, such an object keeps its shape during the simulation which means that the distances between particles do not change.

Let Ω be a rigid object composed by n particles. A local frame is attached to Ω, with R a 3x3 matrix describing its rotation and T a vector describing its translation. The origin of the local frame is given by the center of mass G of Ω. The world position P i of particle i in the object can be expressed in terms of its local frame position P 0 i as P i = RP 0 i + T. Furthermore, the rigid object has a linear velocity V and an angular velocity W. The linear acceleration of the rigid object, V, is computed as follows:

V = n ∑ i=1 F i / n ∑ i=1 m i (1)
where F i is the total of forces applied on the particle i and m i its mass. Also, the angular acceleration of the rigid object, Ẇ, is:

Ẇ = I -1 (τ + (IW) × W) ( 2
)
where I is the inertia matrix associated with the global frame and τ the total torque applied on Ω. τ is computed from the following expression:

τ = n ∑ i=1 GP i × F i (3)
From V and Ẇ, the new linear and angular velocities of Ω are computed by an Euler integration scheme in order to update the new position and orientation of Ω.

Linking rigid and elastic objects

A specific algorithm using direct constraint projection has been developed to simulate elastic or contractile objects fixed to a rigid object, such as a muscle fixed on a bone.

The method constrains the particles located at the body surfaces to keep the same position during the simulation. More precisely, the particles of the interface are cloned and the movements of the elastic particles are constrained by those of the rigid particles. Let Ωe (resp. Ωr) be an elastic (resp. a rigid) body. To handle the interaction between Ωe and Ωr, we associate to each particle Pe ∈ Ωe (resp. Pr ∈ Ωr) a clone. During simulation, the original particles and their clone are kept at the same position. Our collision algorithm will preserve the relative position of Pe and Pr so that the attachment between Ωe and Ωr are locally and transiently preserved. The different steps of the algorithm are the following:

1. compute the forces Fe applied on Pe (the sum of external and internal elastic forces).

2. Calculate the forces Fr applied on Pr (the sum of external forces) and add the forces transmitted by Pe in F i : F i = Fr + Fe. 3. Using (1), (2) and (3) given in section 3.2, compute the new orientation and position of Ωr. Update the position of Pr. 4. Constrain Pe to keep the same position than Pr (i.e. displacement of Pe is forced to be the displacement of Pr).

The algorithm of this method is detailed below for the general case (Alg. 2).

Algorithm 2: linking elastic and rigid objects.

foreach

Sticking and adherence

The proposed algorithm to model a sticky state is given in this section. The first paragraph is devoted to the collision detection between deformable objects. The method to provide a sticky state is then detailed in a second paragraph.

Collision detection

We assume that the object surfaces are represented by triangular meshes. The algorithm proceeds in two stages. First, a coarse detection selects nc pairs of objects potentially in contact. Spheres are classically used as a rough approximation of the objects shape. Second, the particles of these objects involved in the contact are identified. These particles are labelled as Colliding Particles (CPs). The set of actual CPs is determined using a refinement of classical detection algorithm (proximity and interference detection are performed).

Collision response

Our method is based on the creation of a virtual adhesive object in the zone where the bodies are in contact. We consider the adherence as a quasi static state in which a dynamic set c The Eurographics Association 2007.

of CPs are temporarily glued together and act as a rigid body linked to the elastic body in contact (see paragraph 3.2). This guarantees no further interpenetration in the considered zone (due to the projection algorithm explained in the paragraph 3.3). However, unlike [START_REF] Provot | Collision and self-collision handling in cloth model dedicated to design garments[END_REF] or [BFA02] where zones of impact provide non-sliding contacts, our method allows several types of contacts such as rolling, sliding and sticking. Moreover, the link between the objects allows an exchange of forces to take place. Indeed, the elastic forces generated by the deformable objects are transmitted to the virtual adhesive. In return, the resultant of these forces and torques are applied on the CPs. This section is organized as follows: the first paragraph introduces the concept of zone of contact and presents the method to simulate a sticky state. Then, the dynamics of the zones of contact is explained. Finally, the algorithm of the method is detailed.

Zone of contact, clones and virtual adhesive

For each pair of colliding objects, there is one or several zones of contact. A zone of contact is defined by a set of CPs that was detected (Fig. 1a). For each CP of the zone of contact a clone is created. The clone particles are then grouped together. The resulting group, named virtual adhesive (Fig. 1b), is considered as a rigid object which motion provides the collision response.

Clone Each CP belonging to the considered zone of contact is cloned. The clone particle has the same position, mass and velocity as the CP and the two particles are linked as defined in paragraph 3.3. Thus, a CP is constrained to be at the same position as its virtual clone whatever the forces applied on it (Fig. 1b).

Virtual adhesive For each zone of contact a virtual adhesive, composed by virtual clones, is created which allows the simulation of a sticky state between colliding objects in this particular zone (Fig. 1b). In order to freeze the distances between the clone particles, regardless of their origin, the virtual adhesive is modeled as a rigid object. Therefore the contact between objects are maintened without penetration. The virtual adhesive motion depends only on the forces acting on both side of the collision.

Rigid motion During the simulation, all forces applied to the CPs are directly transmitted to their clone. These forces are used to update the virtual adhesive force and torque, which in turn are involved in the rigid object motion (paragraph 3.2 and Fig. 1c, 1d and1e). The colliding objects are temporarily frozen in the considered zone of contact and no further interpenetration is possible (Fig. 1f).

Virtual adhesive population dynamics

As stated before, the interaction between two deformable objects can change the zone of contact and, therefore, a change in the CP population. For example, two colliding objects can part away from one another, decreasing the number of CPs. Alternatively, if the objects agglutinate more and more, the number of CPs will increase. This section describes the method that allows these changes.

When a CP is identified, it has to be assigned to only one virtual adhesive. Three cases can occur. First case, the CP is isolated and a new virtual adhesive has to be created, with the CP being its first element. Second case, the CP has at least one neighbouring particle already in one virtual adhesive. The CP is thus added to this closest virtual adhesive. Third case, there are two or more virtual adhesives close to the CP, i.e. the CP has more than one neighbouring particle in more than one distinct virtual adhesive. All these virtual adhesives are considered to have grown until they became adjacent; they must be merged. In this case, the CP is added to this new larger virtual adhesive.

At each time step, each CP computes the theoretical displacement, Ue, that would have result from elastic internal forces if it were free from any virtual adhesive. This displacement is then compared to the imposed displacement, U VA , resulting from the virtual adhesive motion. If Ue and U VA are in opposite directions, the norm of their difference is compared to a given parameter κ adh . When this norm is higher than κ adh , the clone is destroyed. As a result, the CP becomes free to move again. The value of κ adh controls the degree of adherence. If κ adh is small the adherence is small: only a small difference between Ue and U VA is required to free the particle from the virtual adhesive.

When a virtual adhesive is composed of less than two CPs, it is no more needed and it has to be destroyed (there are no more collisions). Any virtual adhesive can recruit or loose particles according to the relative motions of the different objects in contact. The value of κ adh controls the virtual adhesive population dynamics allowing rolling, sliding or sticking.

General algorithm

The method proceeds in five consecutive stages (Alg. 3 and Fig. 1):

Stage 1 identifies a set of CPs from each pair of selected colliding objects (paragraph 4.1), defining the zone of contact (Fig. 1a). Stage 2 clones each CP. From these clone particles, a virtual adhesive is created for each zone of contact (paragraph 4.2.1 and Fig. 1b). Stage 3 computes the internal forces applied on each particle of the objects. These forces are summed (including elastic forces) and transmitted from the CPs to their clone (Fig. 1c). The position and orientation of the virtual adhesive are then computed by using these transmitted forces as input forces for the rigid body motion algorithm (paragraph 3.2 and Fig. 1d and1e). In marked contrast to approaches that consist of adding energy to handle collisions [Mir96, PMT03, CBP05], our method uses the existing forces to resolve contacts. Indeed, the elastic internal forces applied to the CPs are transmitted to their clone. These forces are then used by the virtual adhesive to compute its new position and orientation by using equations (1), (2) and (3) (rigid body motion). The collision response for all CPs are then given by the new position of their clone. As a result, the contacts are resolved without introducing additional energy in the system that is to say simply using the forces in presence.

Simulations

All simulations described in this section are performed on an Intel Xeon 5140 2.33 GHz and compiled with gcc v4.1 whithout any particular effort to optimize the code. The method is implemented in C++ and the different models are manipulated in PML [START_REF] Chabanas | Physical model language: Towards a unified representation for continuous and discrete models[END_REF]. All the elastic objects are modeled by particles on their external surfaces, organized in triangular meshes, and are satisfying the volume preservation constraint. In a first experiment, six deformable spheres fall onto an elastic plane due to gravity. The bodies deform on impact and end up embedded in the plane. The scene consists of 693 particles and 1280 triangles and the average time for collision detection and response is 549ms. As shown in Fig. 2, our approach provides an interpenetration-free situation between the objects. There are no interpenetrations between the spheres and the bottom sphere does not intersect with the plane. Note that the differences between the shape of the objects are due to elasticity and volume preservation.

An experiment that compares the collision between two types of sphere is shown in Fig. 3. At the top, both spheres are very adhesive while the spheres on the bottom are less sticky. For 84 particles and 160 triangles the average time for collision detection and response is bellow 10ms for the adhesive spheres and bellow 2ms for the others. At the time of impact, the spheres are colliding and slide on one another. Fig. 3 shows the time series of the difference of behavior between the two types of spheres. Note that in the last screenshot the very adhesive objects are still in contact (top right), while the less adhesive bodies are already in a collision-free situation (bottom right). The number of CPs is given for each screenshot. The experiment shows that the method allows the modeling of different sticky states, depending on κ adh .

A deformable ball rolling onto an elastic plane is simulated in Fig. 4 . The scene consists of 217 particles and 368 triangles. The average time for collision detection and response is 33ms. Fig. 4 shows five frames of the dynamic sequence. Note that the deformation of the plane is due to the sticking contact with the ball and the imperfect rolling motion is due to adhesion between the ball and the plane.

In yet another experiment, our approach is tested with more complex objects. The result of the collision between four bunnies is illustrated in Fig. 5. The method still efficiently and correctly resolves all sticky contacts. This test demonstrates that our method can be applied to a wide range of objects and is not limited to spherical bodies. The scene has 567 particles and 1040 triangles, the average time for collision detection and response is 576ms.

In a last test, an experiment is conducted with a complex scene mimicking cell agregation. An illustration of the aggregate behavior of a group of cells is given in Fig. 6. Twenty five incompressible deformable spheres are attracted by the center of the scene, mimicking the source of an attractive chemical signal, little by little colliding and regrouping into a single aggregate. The scene consists of 1050 particles and 2000 triangles. The average time for collision detection and response is 777ms. The percentage of CPs relative to the total number of particles is shown for each step in the bottom slider. The last screenshot in transparency demonstrates that the method still efficiently and correctly resolves contacts: there is no empty space inside the aggregate.

Conclusion and future work

We present a specific method for the collision, detection and response between deformable objects in the adhesive situation. The method simulates several types of contacts including rolling, sliding and sticking. The key idea is to detect the zones of contact between colliding object and to create virtual bodies in these zones. Each virtual adhesive is composed by particles locally in contact, using virtual clones.

We illustrated our approach with several experiments showing that the method efficiently and correctly resolves contacts. This approach can be applied to objects with a wide range of shapes, although our collision detection stage is only adapted for spherical shape interaction. Finally, the method allows the modeling of several sticky states, which make it suitable for biological applications as well as surgical simulators. proach. Also our method needs to be optimized to not specifically target spherical shapes. A problem that has not been addressed yet is the self-collision. Future work also includes the detection and response for such situations. Finally, we will use the method by Spillmann and Teschner [START_REF] Spillmann | Contact surface computation for coarsely sampled deformable objects[END_REF] to improve the management of interpenetrations between bodies.

Acknowledgment

Figure 1 :

 1 Figure 1: Simulation of two colliding deformable objects: (a) zone of contact defined by the list of five Colliding Particles (CPs); (b) virtual adhesive corresponding to the contact zone, its virtual adhesive particles in green are clones of the five CPs; (c) separated elastic internal forces computation; (d) elastic forces transmitted from the CPs to their clones are used to compute the rigid body motion of the virtual adhesive; (e) resulting torque and force applied to the virtual adhesive;(f) movement of the deformable objects and the virtual adhesive due to the forces. The CPs are constrained to be at the same position than their clone.

Figure 2 :

 2 Figure 2: Incompressible deformable spheres falling onto an elastic plane due to the gravity.

Figure 4 :

 4 Figure 4: An elastic ball rolling onto an elastic plane, from top to bottom: t=0s, t=40s, t=80s, t=120s and t=160s. Note the deformation of the plane due to the sticky contact with the ball.

Figure 3 :

 3 Figure3: Collision between two very adhesive objects (top, high κ adh ) and two less sticky objects (bottom, low κ adh ). On each series, the object of the top is attracted downward, while the object of the bottom is attracted updward, resulting in a collision during the simulation. The figure illustrates the difference of configuration between two values of κ adh at four steps of the simulation. For each screenshot, the number of CPs (N) is specified.

Figure 5 :

 5 Figure 5: Collision between four "sticky" bunnies.

  Currently, we are working on the acceleration of the collision detection step based on an adaptative hierarchical apc The Eurographics Association 2007.

Figure 6 :

 6 Figure 6: Simulation of the aggregation of cells: twenty five incompressible deformable spheres move towards an external source of chemical signal. For each screenshot the percentage of CPs according to the total number of particles is shown in the bottom part. From top to bottom and left to right: t=0s, t=100s, t=200s and t=400s. The last frame illustrates in transparency the efficiency of the contact resolution.

A

  .C. acknowledges a fellowship from the French Ministry of Research and the Rhône-Alpes Institute of Complex Systems (IXXI), France. References [BFA02] BRIDSON R., FEDKIW R., ANDERSON J.: Robust treatment of collisions, contact and friction for cloth animation. In SIGGRAPH '02: Proceedings of the 29th annual conference on Computer graphics and interactive techniques (New York, NY, USA, 2002), ACM Press, pp. 594-603. [BJ07] BARBI Č J., JAMES D.: Time-critical distributed contact for 6-dof haptic rendering of adaptively sampled reduced deformable models. In SCA '07: Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation (2007), Eurographics Association, pp. 171-180. [BMB * 07] BERRO J., MICHELOT A., BLANCHOIN L., KOVAR D., MARTIEL J.: Attachment conditions control c The Eurographics Association 2007.

  Algorithm 1: general algorithm of the model.

	foreach object of the model do
	foreach particle of the object do
	applyBoundaryConditions()
	end
	end
	foreach object of the model do
	foreach particle of the object do
	computeForces()
	redistributeInternalForces()
	end
	end
	foreach object of the model do
	foreach particle of the object do
	sumForces()
	computePosition()
	applyConstraints()
	end
	end

c The Eurographics Association 2007.

c The Eurographics Association 2007.


Michel Fliess
email: michel.fliess@polytechnique.edu

Stefan Fuchshumer
email: stefan.fuchshumer@vatron.com

Markus Schöberl
email: markus.schoeberl@jku.at&kurt.schlacher@jku.at

Kurt Schlacher

Hebertt Sira-Ramirez

Hebertt Sira-Ramírez

An introduction to algebraic discrete-time linear parametric identification with a concrete application

Keywords: Discrete-time linear systems, z-transform, bilinear Tustin transform, parameter identification, on-line identifier, rapid sampling, computer algebra. MOTS-CLÉS : Systèmes linéaires discrets, transformation en z, transformation bilinéaire de Tustin, identification paramétrique, identifieur en ligne, échantillonnage rapide, calcul formel

come

An introduction to algebraic discrete-time linear parametric identification with a concrete application

Introduction

This article, the origin of which is a recent conference communication [START_REF] Fliess | Analyse non standard du bruit[END_REF], further develops recent works on discrete-time parametric identification (see [START_REF] Sira-Ramírez | On discrete-time uncertain visual based control of planar manipulators: an on-line algebraic identification approach[END_REF][START_REF] Fliess | On-line discrete-time linear parametric identification: an algebraic approach[END_REF][START_REF] Fuchshumer | Algebraic linear identification, modelling, and applications of flatness-based control[END_REF]. It is a counterpart, augmented with application issues, of [START_REF] Fliess | An algebraic framework for linear identification[END_REF]Fliess et al., 2008b), which permits for linear time-invariant continuous-time systems, thanks to algebraic methods, to achieve, -on-line parametric identification, -robustness with respect to noisy data without knowing the statistical properties of the corrupting noises (see [START_REF] Fliess | Analyse non standard du bruit[END_REF] for further details).

In accordance to the continuous-time framework, the operational representation of the discrete-time constant linear system (in the z-domain) is considered. Initial conditions are allowed to being ignored by taking derivatives with respect to the shift operator z. To determine the unknown system parameters, subsequent iterated summations of the discrete-time counterpart of the resulting operational equation are carried out to set up a system of linear equations, referred to as linear identifier. The presentation is evolved for a general n-th order discrete-time constant linear dynamics. To cope with measurement noise, besides the possibility of a straightforward incorporation of linear filters, the setup of an over-determined system of linear equations by means of additional iterated summations qualifies to be suitable.

On the basis of a fifth-order model of a drive-train, which is available as a laboratory experiment, the problem of inaccurate estimation of those system zeros, which have only minor effect on the system response, and, hence, are difficult to estimate in presence of noise, is illustrated. Additionally, it is found that the identifier parameterized in terms of the z-domain parameters, i.e., the coefficients of the system's difference equation, exhibits increasingly poor numerical condition emerging with decreasing sampling times. This numerical issue of the z-domain setting might become apparent by reflecting the relation z i = exp (s i T a ) between the poles s i of the continuous-time system and the poles z i of the according discrete-time representation. Hence, with decreasing sampling time T a , the poles z i approach the point z = 1. In order to overcome this numerical problem, a suitable re-parametrization of the identifier is sought for and found in terms of the bilinear Tustin transform, also referred to as q-transform for short in linear systems theory (see, e.g., [START_REF] Middleton | Digital Control and Estimation: A Unified Approach[END_REF] for a related transform). Again referring to the drive-train example it is shown that the q-domain setting of the identifier does not experience numerical deficiencies in case of small sampling times.

To cope with the problem of inaccurate estimation of "inessential" zeros, the idea to discarding (or pre-setting) those zeros, based on a priori knowledge from modeling, is proposed. This approach is first discussed on the basis of the q-domain setting and then transferred to the z-domain framework accordingly. The motivation for discussing this idea for the q-domain case first simply is that, by virtue of the close sim-ilarity to the continuous-time domain, things might be particularly apparent from the control engineer's point of view. This idea of pre-setting those non-essential zeros to obtain a "reduced-order" linear identifier is shown to provide particular attractiveness. The corresponding computer programs are available.

The accuracy of our identification results is demonstrated by inspecting the tight correspondence of measurement and simulation of a drive-train control application, with the controller designed on the basis of the identified dynamics using a simple loop-shaping procedure (see, e.g., [START_REF] Gausch | Digitale Regelkreise[END_REF][START_REF] Longchamp | Commande numérique des systèmes dynamiques[END_REF]).

A linear identifier for n-th order discrete-time SISO LTI systems

Consider a linear time-invariant discrete-time SISO system of order n n i=0

a i y k+i = m i=0 b i u k+i , m ≤ n, a n = 1 [1]
For estimating the unknown parameters

λ T = [a 0 , . . . , a n-1 , b 0 , . . . , b m ] [2]
we start with the z-domain representation of system [1], which yields an identification approach motivated by [START_REF] Fliess | An algebraic framework for linear identification[END_REF]Fliess et al., 2008b).

Consider a sequence (f k ), where

f k = 0 for k < 0. The formal Laurent series f z (z) = ∞ i=0 f i z -i , called the z-transform of (f k ), is written f z ¡(f k ) for short. With (f k+i ) z i f z - i-1 j=0 f j z -j , and A (z) = z n + n-1 i=0 a i z i , B (z) = m i=0 b i z i , the z-domain counterpart of Equation [1] is A (z) y z - n i=0 a i z i i-1 j=0 y j z -j = B (z) u z - m i=0 b i z i i-1 j=0 u j z -j [3]
and the associated transfer function is

G (z) = y z (z) u z (z) = B (z) A (z) [4]
Correspondingly to the continuous-time framework of [START_REF] Fliess | An algebraic framework for linear identification[END_REF]Fliess et al., 2008b), derivatives with respect to z, of order n ≥ n + 1, are taken on both sides of Equation [3] in order to eliminate the initial conditions y j , j = 0, . . . n -1, and u j , j = 0, . . . m -1. (Accordingly, one might think, e.g., of first dividing both sides of Equation [3] by z, followed by an n times differentiation with respect to z, in order to meet this objective.) By carrying out the derivative of order n ≥ n + 1 on both sides of Equation [3] (the detailed calculations are given in the appendix) and finally transferring the result back to the discrete-time domain, we end up with

n-1 s=0 (k -n + s) y k + n-1 i=0 y k-n+i a i - m i=0 u k-n+i b i = 0 [5]
For determining the set λ of parameters, Equation [5] is subject to N ≥ n+m fold iterated summation (correspondingly to the iterated integrals of [START_REF] Fliess | An algebraic framework for linear identification[END_REF]) to assemble a set of (N + 1) linear equations, P (k) λ = Q (k), also referred to as an (on-line) linear identifier. So, with χ (k) = n-1 s=0 (k -n + s), the first row of P and Q is given as

P 1,i+1 (k) = χ (k) y k-n+i , i = 0, . . . , n -1, P 1,n+j+1 (k) = -χ (k) u k-n+j , j = 0, . . . , m, Q 1 (k) = -χ (k) y k . The subsequent rows read P 1+l,i (k) =
k-1 r=0 P l,i (r), l = 1, . . . , N , i = 1, . . . , n + m + 1, with an according expression for Q. The reason for optionally setting up an over-determined system of linear equations (by choosing N > n + m) is the following: it provides an improved accuracy for the estimates in presence of noisy signals. Section 5 will deal with some related issues.

Re-parameterizing the problem via application of the Tustin transform

Given a continuous-time system G (s) with poles s i , then the poles of the according discrete-time system G (z) are located at z i = exp (s i T a ). Thus, with a decreasing sampling time T a , the poles z i approach 1. The objective of the following discussions is to introduce a suitable re-parametrization of the proposed identification method, namely via application of the bilinear Tustin transform C → C,

z = 1 + q/Ω 0 1 -q/Ω 0 , q = Ω 0 z -1 z + 1 , Ω 0 = 2 T a [6]
in order to avoid poor numerical conditioning emerging with decreasing sampling times. Section 6.1 will illustrate these (numerical) issues on the basis of a "drivetrain" example.

Let

Ǵ# (q) = G (z)| z= 1+q/Ω 0 1-q/Ω 0 [7]
denote the q-domain transfer function, with the prime arranged to indicate the representation in terms of λ. Clearly, an (n -m)-fold zero at q = Ω 0 occurs. Next, arrange a change of the parameters to re-cast Ǵ# (q) as

G # (q) = (1 -q/Ω 0 ) n-m m i=0 B i q i 1 + n i=1 A i q i [8]
represented in terms of the parameters {A i , B i }. Notice that, given a continuous-time system G (s), the approximation G # (jΩ) ≈ G (jω), with Ω = (2/T a ) tan (ωT a /2), holds for |ωT a | ≪ 1. Let s = a be a pole of G (s), then the according pole of G # (q) is ã = Ω 0 tanh (a/Ω 0 ). This observation draws a close link between the q-domain and the s-domain scenario.

The following idea will be crucial for additionally increasing the performance of the identifier (to be illustrated in Section 6.1). Being particularly apparent in the qdomain setting from the control engineer's point of view, one may discard, except for the (n -m)-fold zero at Ω 0 , those zeros of G # (q) which only have minor influence on the system response, and, hence, are difficult to estimate in presence of noise. This idea of discarding the "inessential" zeros a priori (based on knowledge obtained from modeling) will be accordingly transferred to the z-domain setting in Section 4.

So, the point of departure for the q-domain setting of the identification approach, incorporating the (optional) feature of discarding "inessential" zeros, is set as

G# (q) = (1 -q/Ω 0 ) n-m m i=0 Bi q i 1 + n i=1 A i q i , m ≤ m [9]
with the tilde indicating the approximation of G # (q). Next, we will trace back to the z-domain solution, and re-cast [5] in terms of the q-domain parameters

Λ T = A 1 , . . . , A n , B0 , . . . , B m [10]
To this end, let us first determine the parameter transformation. The z-domain transfer function (with

a n = 1 in general) according to G# (q) reads G (z) = 2 n-m m i=0 Bi Ω i 0 (z -1) i (z + 1) m-i (z + 1) n + n i=1 A i Ω i 0 (z -1) i (z + 1) n-i = m s=0 b s z s n s=0 a s z s [11]
By referring to the binomial theorem, we obtain the relations between the (n + m + 2) coefficients {a s , b s } of [11] and the (n + m + 1) parameters Λ as

b s = m i=0 Π b s,i Bi , s = 0, . . . , m a s = n s + n i=1 Π a s,i A i , s = 0, . . . , n [12] with Π b s,i = 2 n-m (-Ω 0 ) i s j=0 (-1) s+j i s -j m -i j , Π a s,i = (-Ω 0 ) i s j=0 (-1) s+j i s -j n -i j [13]
Plugging the transformation [12, 13] into [5] (notice that [5] has to be slightly adapted so as to allow for arbitrary a n ), we finally find

χ (k) n r=0 n r y k-n+r + n i=1 n r=0 y k-n+r Π a r,i A i - - m i=0 m r=0 u k-n+r Π b r,i Bi = 0 [14]
The procedure of setting up the linear identifier by means of iterated summation of [14], and the solution for Λ as well, proceeds as discussed in Section 2.

The z-domain setting -cont'd

The idea to (optionally) a priori discarding "inessential" zeros will now be equivalently applied to the z-domain setting, i.e., the parametrization in terms of λ, continuing Section 2. Clearly, by [6], discarding certain zeros of G # (q), i.e., shifting those zeros to infinity, corresponds to placing the associated zeros of

G (z) at z = -1. Now, introduce the counterpart of [9], G (z) = (z + 1) m-m m i=0 bi z i n i=0 a i z i = m i=0 b i z i n i=0 a i z i

Robustness with respect to noisy data

In presence of noisy signals, the composition of the identifier as a "square" system of linear equations, i.e., N = n+m or N = n+ m respectively, via iterated summation of [5] or [14], might not yield accurate (or even appropriate) estimates. For instance, this is the case for the drive-train example to be investigated in Section 6.1. We will address here some possibilities (which might of course be combined) for improving the performance of the proposed identification method in presence of noisy signals:

1) instead of taking iterated summations on [5] or [14], i.e., iterated application of 1/ (z -1) in the z-domain, one might more generally iteratively apply (filters of the type) 1/ (z -γ) (or even higher-order ones) for denoising, 2) "invariant filtering" (see, e.g., (Fliess et al., 2008b)): Each row of the linear identifier P (k) λ = Q (k) (or its counterparts of Sections 3 and 4) might be preprocessed with a filter, suitably adjusted utilizing (a priori) knowledge of the system dynamics, 3) instead of setting up a "square" linear identifer, one might choose N ≫ n + m (or N ≫ n + m, resp.) to set up an over-determined system (in the presence of noise) and solve, e.g., for a least-squares solution to (the under-determined linear system) P λ -e = Q, i.e., min λ e 2 2 , which is λ = (P T P ) -1 P T Q, 4) clearly, the "first few" equations of the linear identifier are obviously more affected by noises than the subsequent ones. This is due to the effect of the iterated summations (or the iterated filters of item 1). So, one might think of discarding the "first" α equations from the parameter calculation.

The combination of the items 3 and 4 might read as follows: Set up a number of N = n + m + α + β equations (either by means of iterated summation or filtering), with β denoting the number of additional equations added due to item 3. Hence, from the (1 + N ) equations, the equations no. (α + 1) up to (1 + N ), i.e., a number of 1 + N -α = 1 + n + m + β equations, are used for calculating the (1 + n + m) unknown parameters (in the least-squares sense).

REMARK. -Several aspects of the above discussion are nothing else but a discretetime interpretation of the viewpoint expounded in [START_REF] Fliess | Analyse non standard du bruit[END_REF], where noises are considered as highly fluctuating, or oscillating, phenomena. Let us emphasize once more that we do not need any statistical knowledge of the noises.

REMARK. -See [START_REF] Mboup | New algebraic estimation techniques in signal processing and control[END_REF][START_REF] Mboup | Parameter estimation via differential algebra and operational calculus[END_REF] for a most illuminating comparison between continuous-time algebraic identification methods and least-square techniques.

Application and discussion

Identification

To illustrate the behavior of the presented approach to discrete-time linear systems identification, and, in particular, to reveal some subtleties involved, we will finally discuss a selected application available as a laboratory experiment.

Θ 1 Θ 2 Θ 3 c 1 c 2 ω permanent-magnet dc motor u Figure 1. The lab model "drive-train"
Consider the model of a drive train as depicted in Figure 1. The parameters of the lab setup are as follows: L A = 896 µH (armature inductance), R A = 6.38 Ω (armature resistance), k m = 41 • 10 -3 N m/A (torque constant), c 1 = c 2 = 1.72 • 10 -3 N m/rad (spring coefficients), Θ 1 = 25.65 • 10 -6 kgm 2 , Θ 2 = 6.44 • 10 -6 kgm 2 , Θ 3 = 5.1 • 10 -6 kgm 2 (moments of inertia of the rotors), and d 1 = 3.98 • 10 -6 N ms, d 2 = 0.92 • 10 -6 N ms, d 3 = 2.4 • 10 -6 N ms (coefficients of viscous friction, related to the bearings of the respective rotors). By discarding the dynamics related to the electrical subsystem in the sense of a singular perturbation point of view, the dynamics of the drive-train are obtained as G (s) = ω/û,

G (s) = V 1 + s 8.1 1 + 2ξ1s ω1 + s 2 ω 2 1 1 + 2ξ2s ω2 + s 2 ω 2 2 [16]
with V = 23.7, ω 1 = 12.4, ω 2 = 27.7, ξ 1 = 0.1 and ξ 2 = 0.0083. Thus n = 5, m = 4 for the according discrete-time system G (z).

The following examples (with n = n + 1) provide different case studies regarding the choice of the parameter m ≤ m (for both the qand z-domain setting) and the sampling time T a . In order to cope with noise, the combination of the items 3 and 4 of Section 5 is applied, with α = 7 and β = 10. Numerous case studies showed that for this example, the use of the iterated summations (i.e., γ = 1) is clearly advantageous compared to the choice γ = 1 of item 1 (of Section 5). The "invariant filtering approach" of item 2 was found not to give significant further improvements (to the setting as introduced above), so it is not applied within the following case studies.

All equations of the identifier are normalized (by dividing by the maximum absolute entry of P of the respective row) to improve the numerical conditioning. The linear on-line identifiers start at t = 0, and the pole-zero plots and Bode diagrams of the identified zand q-domain transfer functions given in the figures are due to λ and Λ, evaluated at the final time t end = 1.35s. The pole-zero plots of the nominal dynamics G (z), or G # (q) resp., associated to G (s) as given above, will always be displayed in blue color.

Example 1 Set m = m, T a = 10 ms. The simulation results given in the Figures 2 and3 are associated with the observation that, in presence of noise (chosen as colored), the estimation of the zeros which have minor influence on the system response (see the subplot containing the system output ω) is very poor. Additionally, it is found that, emerging with decreasing sampling times, the numerical conditioning of the linear identifier parameterized in terms of the z-domain parameters λ becomes increasingly worse, in contrast to the q-domain setting. To illustrate this, the conditioning numbers of P (i.e., the ratio of the largest singular value of P to the smallest) for the zand q-domain identifier, evaluated at the final time t end , for different sampling times, are given in Table 1 (cols. 2 and3).

In order to show that these conditioning numbers are only weakly affected by the noise added to the output signal, the conditioning numbers obtained by noise-free simulations are also given in Table 1 (cols. 4 and 5).

T a [ms]

z q z q 10 3.0e10 1.7e10 3.0e10 1.7e10 5 5.1e11 8.0e9 5.2e11 7.8e9 2 2.9e13 1.2e10 2.9e13 1.2e10 1 5.4e14 1.4e10 5.4e14 1.4e10 0.5 9.0e15 1.5e10 8.9e15 1.5e10

Table 1. (see Example 1) Conditioning numbers of P

The first observation of Example 1 (inappropriate estimation of the zeros in presence of noise), associated with the knowledge obtained from modeling, is the motivation for applying the approximations of the numerators as proposed in Sections 3, 4. Notice that, by inspecting for instance the Bode plots of the nominal drive-train dynamics in Figure 3, those mentioned zeros only affect the frequency response in the frequency domain with the magnitude G # (jΩ) located significantly below the zero dB line. The idea of discarding those zeros is addressed in the following example.

Example 2 Set m = 0 and T a = 10 ms for the simulation results given in the Figures 4 and 5. The approach for discarding those zeros having negligible effect on the system response is seen to be appropriate (see in particular the subplot containing the system output ω). The estimation of the poles is again accurate. Additionally, due to the reduced number of parameters, this approach has also the advantage of having better numerical conditioning compared to the case m = m of the previous example, illustrated by the condition numbers given in Table 2. Analogously to Table 1, cols. 4

T a [ms] z q z q 10 1.8e7 2.3e7 1.9e7 2.3e7 5 6.0e8 2.7e7 5.8e8 2.7e7 2 3.5e10 3.0e7 3.5e10 3.0e7 1 6.4e11 3.1e7 6.4e11 3.1e7 0.5 1.1e13 3.2e7 1.1e13 3.2e7

Ω = Ω 0 tan (ω/Ω 0 ) [rad/s]
q-ident. z-ident. nom.

Figure 3. (see Example 1) Figure 2 cont'd: Bode plots

Let us end these case studies with some considerations on the counterparts of the Examples 1 and 2 by invoking a standard least-squares (LS) method in the zand the q-domain. We also present a comparison with an instrumental-variable (IV4) identification method (see, e.g., [START_REF] Ljung | System Identification -Theory for the User[END_REF][START_REF] Söderström | System Identification[END_REF]) in the z-domain, since by assumption the equation error is not white noise.

Example 4 (standard least-squares (LS) method; m = m, T a = 10 ms). The results obtained with the same measurement data as given in Figure 6 are displayed in Figure 8. They indicate poor performance (getting worse with decreasing sampling times), even for the estimation of the pole-pair related to the "slow eigenfrequency".

Example 5 (standard LS method; m = 0, T a = 10 ms). This setting is found not to give any clear improvements compared to the choice m = m of Example 4, thus, the associated graphics are not displayed.

REMARK. -In order to achieve suitable results with the LS identification (in the case of the drive-train example), it is advisable to increase the sampling time, say, e.g., to T a = 50 ms, additionally to articulately increasing the "observation time span".

-2 REMARK. -The use of colored noises is a further confirmation that our denoising techniques are not limited to classic Gaussian white noises 1 .

-1 0 1 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 u [V] t [s] -120 -80 -40 0 40 80 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ω [rad/s] t [s] ω sim. q-ident. z-ident.
Ω = Ω 0 tan (ω/Ω 0 ) [rad/s] q-ident. z-ident. nom.

Control application

Based on the identification result of the drive-train lab model displayed in the Figures 6 and 7, a controller R (for the one-degree-of-freedom standard control scheme) is now designed by means of a loop-shaping procedure, i.e., a simple graphical method based on Bode plots (see, e.g., [START_REF] Gausch | Digitale Regelkreise[END_REF][START_REF] Longchamp | Commande numérique des systèmes dynamiques[END_REF]). The design is carried out in the q-domain. To this end, the transfer function (i.e., the amplitude 1. See [START_REF] Fliess | Analyse non standard du bruit[END_REF] and the references therein for confirmation in various areas including signal processing (see also [START_REF] Neves | Algebraic parameter estimation of damped exponentials[END_REF][START_REF] Trapero | An algebraic frequency estimator for a biased and noisy sinusoidal signal[END_REF][START_REF] Trapero | On the algebraic identification of the frequencies, amplitudes and phases of two sinusoidal signals from their noisy sums[END_REF]). The following control design may be retraced via the Bode plots of Figure 10. Let G # (q) denote the identified q-domain transfer function of the drive-train due to Figure 7. As the first design step, in order to cope with the fairly-damped torsional oscillations of the drive-train system, the according complex-conjugate polepairs -1.16 ± 12.5 √ -1 and -0.76 ± 26.34 √ -1 of G # (q) are compensated for. Let us represent the polynomials associated with those pole-pairs as 1+2ξ i q/ω i +(q/ω i )

2 , i = 1, 2, with (ω 1 = 12.55 rad/s, ξ 1 = 0.0923) for the "slow" eigenfrequency and (ω 2 = 26.35 rad/s, ξ 2 = 0.0289) for the "fast" eigenmode, respectively. Typically, in view of robustness issues regarding performance, it is advisable not to exactly cancel out such pole-pairs 1 + 2ξ i q/ω i + (q/ω i )

2 , but to install (approximate) compensators, as, e.g., The representations in terms of λ and Λ are entitled as "LS-z" and "LS-q" for short with ξi ≥ ξ i ("approximate" compensation as mentioned above takes place for ξi > ξ i ) and 0 < ζ i ≤ 1. Compensators of this type are usually also referred to as Notch filters. In order to illustrate the quality of the identification result, however, we will exactly cancel out the pole-pair related to the "slow" eigenfrequency via N # 1 (q), i.e., set ξ1 = ξ 1 , see also Figure 10. The other pole-pair is compensated approximately by setting ξ2 = 1.5ξ 2 . For both Notch filters of the drive train controller we choose ζ i = 1, i = 1, 2. Finally, as the third part of the controller, in order to achieve steady-state accuracy, we add a transfer function H of PI-type, H # (q) = 0.08 (1 + q/10) /q, to obtain the controller R = HN 1 N 2 . The Bode plot of the open-loop transfer function L = RG is also displayed in Figure 10. The BIBO stability of the closed loop T = L/ (1 + L) is immediately deduced by inspecting the phase angle of L # (jΩ) at the cut-off frequency Ω c = 1.78rad/s, i.e., arg L # (jΩ c ) = -120.7 deg. Measurement Of course, more advanced design techniques, e.g. for decoupling the disturbance and tracking design, might be applied. The primary objective of this simple control example however was to emphasize, by means of measurement and simulation results, the accuracy of the obtained identification result which the control design was based upon.

N # i (q) = 1 + 2 ξi q/ω i + (q/ω i ) 2 1 + 2ζ i q/ω i + (q/ω i ) 2 [17] -2 -1 0 1 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 u [V] t [s]

Consequences of the algebraic setting

The setup of a linear identifier for discrete-time LTI SISO systems, evolving from an algebraic point of view, has been discussed and evaluated on the basis of case studies referring to a fifth-order model of a drive-train. Two different parameterizations have been investigated with regard to their numerical conditioning, with the q-domain setting found to exhibiting significant advantages, which become particularly apparent with decreasing sampling times. This effect is observed independently from whether or not the idea of pre-setting certain zeros, which turned out as useful, is applied.

Though this algebraic approach provides promising results, it is worth mentioning that, clearly, the linear identifier cannot incorporate a priori knowledge on stability. More concretely, for the drive-train example, the pole-pair related to the "fast eigenfrequency" is located very close to the stability margin, hence, in presence of noisy signals, the estimated pole-pair might be found to shift beyond the stability margin.

Finally, it should be mentioned that the discussed method provides easy-toimplement on-line identifiers. A computer-algebra implementation of this approach, associated with notes on implementation issues, is available at http://regpro.mechatronik.uni-linz.ac.at

Improvement of the mathematical formalism

Future publications will associate z-transforms and module theory2 in order to give a more intrinsic algebraic picture of those parametric identification methods and of their extension to multivariable systems. It will thus provide a better understanding of their connections with flatness-based predictive control for discrete-time linear systems ( [START_REF] Fliess | Une approche intrinsèque de la commande prédictive linéaire discrète[END_REF][START_REF] Sira-Ramírez | Differentially Flat Systems[END_REF]) and with structural properties as derived from the module-theoretic standpoint (see, e.g., [START_REF] Bourlès | Systèmes linéaires : de la modélisation à la commande[END_REF][START_REF] Fliess | Une approche intrinsèque de la commande prédictive linéaire discrète[END_REF] and the references therein).

Nonlinear extension

A nonlinear extension of our approach is possible via a discrete-time analogue of (Fliess et al., 2008a).

  indicating the approximation of G (z). The linear map relating the parameters bT = [ b0 . . . b m] to the numerator coefficients b T = [b 0 . . . b m ] is represented as b = Ξ b b. Let λT = [a 0 , . . . , a n-1 , b0 , . . . , b m], then λ = Ξ λ with Ξ = diag (E n×n , Ξ b ), and the linear identifier reads P λ = Q, with P = P Ξ.

Figure 2 .

 2 Figure 2. (see Example 1) Simulation results, m = m, T a = 10 ms. Both zand q-domain setting provide good results (regarding the system dynamics in the "interesting" frequency domain, see the system response and the Bode diagrams of Figure3. However, the estimation of the zeros of both parameterizations is very inappropriate indeed

Figure 4 .

 4 Figure 4. (see Example 2) Simulation results, m = 0, T a = 10 ms. The usefulness of the idea of approximating the numerator in the sense of Sections 3 and 4 is verified, see also Figure 5 for the Bode plots

Figure 5 .

 5 Figure 5. (see Example 2) Figure 4 cont'd: Bode plots

Figure 6 .

 6 Figure 6. (see Example 3) Measurement results, m = 0, T a = 10 ms. See also Figure 7 for the Bode plots

Figure 7 .

 7 Figure 7. (see Example 3) Figure 6 cont'd: Bode plots

Figure 8 .

 8 Figure 8. (see Example 4) Results obtained from measurements by invoking the standard least-squares (LS) method with m = m, T a = 10 ms. The representations in terms of λ and Λ are entitled as "LS-z" and "LS-q" for short

Figure 9 .

 9 Figure 9. (see Example 6) Results obtained from measurements by invoking the standard instrumental variable method (IV4) implemented in MATLAB with m = m, T a = 10 ms

Figure 10 .

 10 Figure 10. Control application: Bode plots of the identification result G # (q) due to the measurements of the lab model (see Figure 7), and the plots of the Notch filters N 1 , N 2 and the open-loop transfer function L = RG

Figure 11 .

 11 Figure 11. Control application: measurement and simulation results of the control system subject to a step change of the reference angular velocity

Table 2 .

 2 (see Example 2) Conditioning numbers of P and 5 display the results obtained with noise-free simulations. Again, the numerical conditioning of the z-domain setting suffers with decreasing sampling times, whereas the parametrization in terms of Λ does not show such problems.

	Example 3 (measurement results) Set m = 0, T a = 10 ms. Figures 6 and 7 depict
	the identification results obtained from measurements of the "drive-train" lab model.
	The angular velocity ω is measured via an incremental encoder.

Remember that operational calculus and module theory were both already employed in[START_REF] Fliess | An algebraic framework for linear identification[END_REF] Fliess et al., 2008b).

Appendix: the detailed calculations of Section 2

The following calculations are given to retrace the appearance of [5] by taking n ≥ n + 1 derivatives on [3] w.r.t. z. Let (•)

(j) = (∂/∂z) j (•), and notice that

Then, taking the derivatives on both sides of [3] yields (Ay z ) n) , and, by virtue of Leibniz' product rule,

noticing that A (j) (z) = 0, j > n, and

To further proceed with [20], and, in particular, in view of facilitating the transformation back to the discrete-time domain, let

Multiply both sides of Equation [20] by z n-n . Then, by re-arrangement of the sums for collecting the parameters a i , i = 0, . . . , n -1 (notice that a n = 1 by [1]), and b i , i = 0, . . . , m, we obtain

Before proceeding with [21], first notice that the identity

holds (with n involved in the left-hand side cancelling out). Then, re-sorting to the expressions of [21] associated to the parameters {a i , b i }, we have z i-n n!