domain: spi chim phys

text:
Grégoire Martouzet

Loren Jørgensen

Yoann Pelet

Anne-Laure Biance

Catherine Barentin
email: catherine.barentin@univ-lyon1.fr

Dynamic arrest during the spreading of a yield stress fluid drop

When a liquid drop is gently deposited on a wetting solid surface, it spreads due to capillary forces until it reaches a thermodynamical equilibrium set by the relative surface energies of the system.

We investigate here experimentally the spreading ability of drops made of yield stress fluids, which flow only if the applied stress is above a finite value. We observe that in this case, after a spreading phase, the motion stops and a well-defined contact angle can be measured. This contact angle depends on the rheological properties of the fluid and in particular on its yield stress, on the drop radius and on the hydrodynamic boundary condition at the surface. These results are quantitatively compared to an analysis showing that, due to the yield stress of the fluid, a mechanical equilibrium is indeed reached which does not correspond to the thermodynamical equilibrium.

I. INTRODUCTION

When a drop of a Newtonian liquid is gently deposited on a smooth substrate, it can spread totally or partially, depending on the nature of the materials (the solid surface, the surrounding gas, the spreading liquid). During the spreading, a so-called contact line is established at the limit of the three phases (here gas, liquid and substrate) and an equilibrium contact angle θ e is reached. θ e satisfies both the minimization of the energy of the system and the force balance at the contact line [1]. In the case of a partial spreading on an ideal surface, θ e is given by Young-Dupré's law: cos θ e = (Γ sg -Γ sl )/Γ, where Γ sg , Γ sl and Γ are the interfacial tensions of the solid/gas, solid/liquid and liquid/gas interfaces, respectively. A dynamic regime precedes equilibrium, during which the contact angle θ decreases and the contact radius R of the drop on the surface increases with time. These dynamics have been extensively studied [START_REF] Tanner | A theory of die-swell[END_REF][START_REF] Voinov | Hydrodynamics of wetting[END_REF][START_REF] Gennes | Wetting: statics and dynamics[END_REF][START_REF] Biance | First steps in the spreading of a liquid droplet[END_REF][6], and the main difficulty has been to model the singularity at the moving contact line, where viscous dissipation diverges. The introduction of a cut-off length and the small-angle and small capillary number approximation lead to the well-established Cox-Voinov law, first derived by Tanner [START_REF] Tanner | The spreading of silicone oil drops on horizontal surfaces[END_REF], that links the contact line velocity with the dynamic contact angle. However, this commonly accepted description is put into question in the case of complex surfaces [START_REF] Lafuma | Superhydrophobic states[END_REF][START_REF] Perrin | Defects at the nanoscale impact contact line motion at all scales[END_REF][START_REF] Andreotti | Statics and dynamics of soft wetting[END_REF], or complex fluids as investigated here.

Different types of complex fluids can be considered. On the one hand, some non-Newtonian fluids flow like liquids whose viscosity depends on the applied shear rate: a decreasing (resp. increasing) viscosity with the shear rate corresponds to a shear-thinning (resp. shear-thickening) behavior. Spreading dynamics of such shear-thinning systems has been studied theoretically [START_REF] Starov | Spreading of non-newtonian liquids over solid substrates[END_REF] and experimentally [START_REF] Rafaï | Spreading of non-newtonian fluids and surfactant solutions on solid surfaces[END_REF] and deviates slightly from what is observed for Newtonian liquids. The rheology of the fluid modifies the spreading dynamics, as expected, but not the equilibrium state, which still satisfies surface energy minimization and Young-Dupré's law. On the other hand, fluids with more complex mechanical responses such as foams, concentrated emulsions, polymer microgels, colloidal suspensions and wet granular materials can be encountered. These materials, called yield-stress fluid (YSF), behave like an elastic solid below a critical stress referred to as the yield stress, and above it, they flow like a liquid. This complex behavior that is both liquid and solid is important for many applications in various industries (food, cosmetics, building) [START_REF] Mezzenga | Understanding foods as soft materials[END_REF][START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF]. More specifically, many processes such as 3D printing, coating or imbibition involve YSFs in contact with surfaces and thus motivate the study of their wetting properties.

In the last decade, studies have examined interfacial properties of YSFs in some situations such as drop formation [START_REF] German | Formation of viscoplastic drops by capillary breakup[END_REF], drop impact [START_REF] Luu | Drop impact of yield-stress fluids[END_REF][START_REF] German | Impact of shear-thinning and yield-stress drops on solid substrates[END_REF][START_REF] German | The spreading behaviour of capillary driven yield-stress drops[END_REF], dip-coating [START_REF] Boujlel | Measuring the surface tension of yield stress fluids[END_REF][START_REF] Maillard | Dip-coating of yield stress fluids[END_REF][START_REF] Gans | Dip-coating of suspensions[END_REF][START_REF] Smit | Stress field inside the bath determines dip coating with yield-stress fluids in cylindrical geometry[END_REF], capillary adhesion [START_REF] Barral | Adhesion of yield stress fluids[END_REF][START_REF] Jørgensen | Yield stress and elasticity influence on surface tension measurements[END_REF] and capillary rise [START_REF] Géraud | Capillary rise of yield-stress fluids[END_REF]. They highlighted the role of the non-relaxed internal stresses in YSF, leading for example to deviations from the classical Jurin's law in capillary rise [START_REF] Géraud | Capillary rise of yield-stress fluids[END_REF] or making the measurement of their surface tension difficult [START_REF] Boujlel | Measuring the surface tension of yield stress fluids[END_REF][START_REF] Jørgensen | Yield stress and elasticity influence on surface tension measurements[END_REF]. But none have considered the spreading of a YSF drop on a solid surface despite the importance of this case both conceptually and practically.

Here we investigate the spreading of a drop made from a YSF on a solid surface and we focus on the influence of the yield stress on the final state, i.e., on the shape of the final drop and on the final contact angle. In particular we address the following question: Does the classical Young-Dupré's law still hold? If not, is it possible to predict it? First, we describe our experimental procedures for monitoring the spreading of drops made of YSF. Then we present our experimental results on the final contact angle observed as a function of the yield stress, which we compare quantitatively to analytical calculations. Finally, a phase diagram showing the different outcomes on the interactions of a YSF drop with a surface is proposed.

II. EXPERIMENTS AND RESULTS

Drops made of various fluids are formed at the edge of a pipette tip. Carbopol gels (ETD2050 and U10 from Lubrizol ), which are dense suspensions of microgels made of cross-linked polyacrylic acid in water, have been used as model YSF. Different concentrations have been considered to achieve different rheological properties, as reported in Tab. I. These rheological properties (yield stress σ y and elastic modulus G ) are measured with a rheometer (Anton Paar Physica MCR 301) in a parallel plate geometry (see Supplemental Material at URL for details on rheological measurements). Yield stresses vary between 0 and 75 Pa. The fluid/gas interfacial tension Γ is measured using a home-made set-up [START_REF] Jørgensen | Yield stress and elasticity influence on surface tension measurements[END_REF], and is around 60 mN/m, for all the fluids considered. Suspensions of non-crosslinked linear polyacrylic acid (PAA) (from PolySciences), which has no yield stress but similar surfaces properties to Carbopol, have been used for comparison. In the following, Carbopol gels will refer to YSF and PAA suspension will refer to fluid without yield stress.

Drops of different initial radii (R 0 = 0.5-2 mm) are gently deposited, with a velocity below 1 mm/s, on a hydrophilic microscope glass slide (from VWR) that is cleaned with an ozone plasma (see Supplemental Material at URL for cleaning procedure). To get smooth substrates, the microscope glass slides are used as received, while for rough ones, the slides are sandblasted. The measured roughness of the substrate is ∼ 20µm, always larger than the microstructure of the YSF used, ensuring a no-slip boundary condition at the solid interface [START_REF] Péméja | Wall slip regimes in jammed suspensions of soft microgels[END_REF]. As shown in Fig. 1a), the experiments are performed in a sealed box saturated with water vapor to avoid evaporation of the drop. The sealed box is placed in a clean room environment to avoid dust deposits on the microscope glass slides during their transfer from the plasma cleaner to the box.

When the drop touches the substrate, it initially spreads. After a few minutes, the motion stops, which indicates that a final state is reached. A side-view picture is taken with a camera (Nikon d700 with a Navitar Zoom 6000 Lens).

As shown in Fig. 1c), different non-regular drop profiles are observed, depending mainly on the yield-stress value and on the drop radius. From these profiles, we extract the final radius (R f ) and final contact angles at the right and left edges of the drop (θ r , θ l ), which are both defined as the greatest slope of the profile ξ(r) on the right (respectively left) side of the drop (see Supplemental Material at URL for the characterisation of contact angles). This extraction is made with a home-made Python code, very similar to ImageJ's Drop Shape Analysis plugin [START_REF] Stalder | A snake-based approach to accurate determination of both contact points and contact angles[END_REF]. The accuracy on the angle is 5 o and 1% for the radius. For symmetric shapes as the ones shown in Fig. 1c) (drops I and II), the final contact angle (θ f ) is defined as the average over θ l and θ r . When the observed shape is clearly asymmetric, characterized by a difference of left and right angles larger than 10%, such as the one shown in Fig. 1c) (drop IV), the final contact angle θ f is not measured and not compared to the model, which assumes a cylindrical symmetry. a * corresponds to values of surface tension that have not been measured in [START_REF] Jørgensen | Yield stress and elasticity influence on surface tension measurements[END_REF] (see Supplemental Material at URL for more details).

Variations of the final contact angles θ f as a function of the final radii R f are reported in Fig. 2a) for the various complex fluids used, and for smooth and rough substrates. On the contrary to what is observed with Newtonian fluids and to what is expected from the classical thermodynamic description, the final contact angle obtained with the YSFs (Carbopol ETD and U10) depends on the drop size; it increases with R f . As a comparison, we also plot θ f as a function of R f for a PAA suspension, which has no yield stress (see green circles, same figure). As expected for such a system, the final contact angle is almost constant with respect to the drop radius. It is also almost constant with the roughness of the substrate. The second important result is that the final contact angle measured with Carbopol (U10) depends on the substrate roughness (see red squares in Fig. 2a)): rough substrates result in contact angles θ f larger of at least 20 • for this YSF. Such dependence of the contact angle with the substrate roughness observed with YSF and more importantly the fact that final contact angles are larger on rough substrates than on smooth ones is opposited to the behavior of Newtonian fluids. Indeed for θ e < 90 • , the contact angle measured with Newtonian fluids on a rough substrate is smaller than the one measured on a smooth substrate due to the gain in surface energy, as described by Wenzel's law [START_REF] Wenzel | Resistance of solid surfaces to wetting by water[END_REF]. Finally, the last important result is the increase of the final contact angle with the yield stress of the fluids, as reported in Fig. 2b), for a given drop size spreading on a rough substrate. σ the applied shear stress, γ the shear rate, K the consistency and n the so-called HB exponent. We then describe the flow profile in the spreading drop, the drop being characterized by its radial height ξ(r). We note R the contact radius of the drop and V the velocity of the contact line, V = dR/dt. Notations are defined in Fig. 1b). Following the derivations for Newtonian fluids spreading on an hydrophilic surface [START_REF] Tanner | A theory of die-swell[END_REF], we assume that (i) the drops are sufficiently small to neglect gravity, (ii) the dynamical contact angle is small so that the lubrication approximation is satisfied and (iii) inertia is negligible so that the drop has a quasi-static equilibrium shape. The last two assumptions are relevant as soon as the dynamic contact angle is smaller than 30 • [2] and the time after the contact is larger than a few milliseconds [START_REF] Biance | First steps in the spreading of a liquid droplet[END_REF]6], respectively.

To determine the spreading dynamics, we first write a momentum balance equation that links the velocity profile within the drop and the pressure gradient. Then, by volume conservation, the pressure gradient is linked to the velocity of the contact line, V . Finally, balancing the work done by the capillary force with the viscous dissipation allows us to close the problem.

More precisely, momentum balance links the shear stress σ at a position (r, z) within the drop to the radial pressure gradient A(r) as: σ(r, z) = A(r) • (ξ(r) -z). Here we assume a free interface. Using the HB law, the relation between the shear rate and the radial velocity γ = ∂v r /∂z and assuming a no slip boundary condition, we calculate the velocity profile within the drop that reads:

           v r (z) = n n + 1 A K 1 n z 1 n +1 y 1 -1 - z z y 1 n +1 for 0 ≤ z ≤ z y v r (z) = n n + 1 A K 1 n z 1 n +1 y for z ≥ z y (1)
Two distinct zones are identified: a sheared zone near the substrate up to an altitude of z y (r) = ξ(r) -σ y /A(r), and a non-sheared zone above that. If ξ(r)A(r) ≤ σ y then z y (r) = 0 and the drop stops.

Using Eq. 1, we then apply the flow rate conservation ξ 0 2πrv(r, z)dz = Q(r) and numerically determine the pressure gradient A(r), assuming a spherical cap for the shape of the drop. Finally, a balance between the capillary and viscous forces allows one to predict the full spreading dynamics. Specifically, the power of the driving capillary force reads P cap = 2πRΓ(cos θ e -cos θ)V whereas the viscous dissipated power in the drop is:

P diss = 2π R 0 ξ a rσ γ dz dr ( 2
)
where a is a cut-off length introduced to regularize the contact line singularity. For Carbopol gels, which are a dense suspension of polymer microgels, we can assume that this cut-off length a is set by the microstructure, i.e. the size of a microgel (∼ 1µm). Note that both powers (P diss , P cap ) depend implicitly on time via the geometrical parameters (R, θ, ξ) which are also time-dependent. Following the methodology of Tanner [START_REF] Tanner | The spreading of silicone oil drops on horizontal surfaces[END_REF], we assume that the dominant part of the viscous dissipation takes place in the corner near the contact line, in r = R, simplifying the problem as a 2D one on a contour length 2πR (see appendix A for details). Taking this assumption, the resulting balance

P cap (V, t) = P diss (V, t
) is solved numerically at each time step that allows to determine V (t). The spreading is then assumed to stop when this balance has no solution any more or when z y = a. The first condition corresponds to vanishing spreading velocity and the second corresponds to the situation where the sheared zone has reached the size of a polymer microgel. From this, we can then determine θ f and R f as a function of the input parameters, such as θ e , Γ, σ y , and the initial volume of the droplet. Results of cos θ e -cos θ f as a function of R f σ y /Γ are reported in Fig. 3.

B. End of the spreading

To get more physical insights, we consider the asymptotic limit near the final state. When the motion is about to stop, the stress contributing to the dissipation, i.e. in the part of the liquid that is flowing, is dominated by the yield stress, so that the stress σ in Eq. 2 can be approximated by σ y . Moreover, in the corner region considered here, we can show that ξ a γ dz = V (see Appendix A). This implies that the viscous dissipation can be approximated by

P diss ∼ σ y V πR 2 .
When balanced with the power done by the capillary forces near the final state, i.e. θ = θ f , the final contact angle can be derived as:

cos θ e -cos θ f ∼ σ y R f Γ . ( 3
)
The characteristic dimensionless number Bc = σ y R f /Γ that compares the yield stress and the capillary pressure, known as the Bingham capillary number [START_REF] Bertola | Wicking with a yield stress fluid[END_REF], is then introduced. Eq. 3 predicts that the larger the yield stress and the drop radius, the larger the final contact angle, as observed experimentally in Fig. 2. In the limit of large Bc, the final contact angle θ f is not defined and the approximations of the model should be revisited, as discussed later.

To test the prediction of the scaling law (Eq. 3), we plot (cos θ e -cos θ f ) obtained experimentally and numerically as a function of Bc = σ y R f /Γ, in Fig. 3. The value of the equilibrium contact angle θ e is here set to 10 • which is the angle measured with the samples exhibiting no yield stress. For small Bingham capillary number (Bc ≤ 1), we find a very good agreement between the scaling law (Eq. 3), the experimental data measured on rough surfaces, and the numerical calculations assuming a no-slip boundary condition. A numerical prefactor of the order of one (0.50 ± 0.03) is obtained when fitting the data with Eq. 3 (slope of the black dotted line in Fig. 3). We also see that the final contact angles measured with Carbopol gel (ETD) on smooth surfaces (open blue diamonds) are well described by the numerical resolution and the asymptotic limit (Eq. 3) suggesting that a no-slip boundary condition is also satisfied for these experiments. This is not the case for the measurements performed with Carbopol gel (U10); the final contact angles measured on smooth surfaces (open red squares) are much smaller than the ones measured with rough surfaces (plain red squares). We attribute this behavior to the existence of a slip at the wall, observed in different situations but with similar (YSF, smooth substrate) couples [START_REF] Péméja | Wall slip regimes in jammed suspensions of soft microgels[END_REF]. In order to take this new boundary condition into account, we introduce in our model a slip velocity at the solid wall, denoted V s . Here, we assume a partial slip boundary condition, i.e., V s < V . Friction laws reported with YSFs [START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF][START_REF] Péméja | Wall slip regimes in jammed suspensions of soft microgels[END_REF][START_REF] Meeker | Slip and flow in soft particle pastes[END_REF][START_REF] Zhang | Wall slip of soft-jammed systems: A generic simple shear process[END_REF] show that the slip velocity V s essentially depends on the stress at the wall σ wall = σ(r, z = a), which a priori depends on the radial position r. However, at the end of the spreading, the stress at the wall tends toward the yield stress σ y so that, in the corner, V -V s can be considered constant in Eq. 2. The dissipative power then scales as P diss ∼ πR 2 σ y (V -V s ) and the balance of powers leads to the following relation:

cos θ e -cos θ f ∼ 1 - V s V σ y R f Γ . (4)
where V s /V corresponds to the the fraction of slippage in the limit of small velocities. Using this new relation (Eq. 4), a very good agreement is found between the measurements made on smooth surfaces (red open squares in Fig. 3) and the scaling law assuming V s /V = 0.66 (red dashed line in Fig. 3). Finally the ability of YSFs to slip depends on the microscopic properties of the couple (YSF, substrate), for example on the microgel size and elasticity and on the substrate roughness. In particular, it has been shown [START_REF] Péméja | Wall slip regimes in jammed suspensions of soft microgels[END_REF] that for a given wall stress and fluid yield stress, Carbopol gel (U10) exhibits slip velocities on smooth surfaces much larger than Carbopol gel (ETD). This explains the difference of spreading behavior measured on smooth surfaces.

IV. DISCUSSION

A crucial outcome from this analysis is that the final contact angle is determined by both the Bingham capillary number Bc and the hydrodynamic boundary conditions at the solid surface. Our analysis is however valid only under three main assumptions we want to discuss here (i) the drop has the shape of a spherical cap (ii) viscous dissipation takes place predominantly near the contact line (iii) gravity is negligible.

The first main assumption of our analysis is that the shape of the drop is a spherical cap. To test the validity range of this assumption, we fit the drop shapes with a spherical cap of a same volume, as shown in Fig. 1c). To quantify the deviation from the spherical shape, we define the quantity δ as:

δ = 1 R f π 0 (ρ sc (α) -ρ(α)) 2 dα, (5)
taking polar coordinates in the picture plane (angle α, radius ρ and the origin is taken in the center of the drop in the plane of the substrate). ρ and ρ sc correspond to the drop and spherical cap profiles respectively. We then report in One can observe that at large final radius R f or large yield stress σ y , δ is significant so that the deviations from the spherical shape cannot be ignored. A frontier between the spherical and non-spherical regimes can be determined, by comparing the Laplace pressure inside the drop which tends to make it spherical to the yield stress which tends to freeze its shape. This balance reads Bc = 1 which is plotted in Fig. 4 (see dashed line). This curve indeed properly delimits the two regions of spherical and non-spherical drops observed experimentally and δ correctly quantifies the deviation from sphericity. Our full analysis is then valid provided that the Bingham capillary number satisfies Bc ≤ 1.

Outside this limit, the history of drop formation and shape has to be considered. In case of conserved sphericity, note that this limit Bc ≤ 1 corresponds also to regimes where the final contact angle is larger than 1, and then where lubrication approximation fails.

The second main assumption is that viscous dissipation is mainly located near the contact line. Similarly to what is observed for Newtonian fluids, we can indeed estimate that near the contact line, the dissipated power integrated over the thickness of the droplet ξ is diverging. In our complex situation, the liquid is sheared only on a height z y ξ, and the dissipated power p per unit length at the contact line then reads, combining Eq. 2 and the complex fluid constitutive law: p ∼ z y K(V /z y ) n+1 . Near the contact line, the height of the droplet tends to zero, so is the sheared zone z y , and consequently the dissipation tends to infinity. As classically done for Newtonian fluids, this non-physical divergence is regularized by introducing the microscopic cut-off length a as detailed before. This assumption is then always valid.

The last assumption of our analysis is that gravity is neglected with respect to surface tension. This has been more

classically analyzed [1] and is governed by the so-called Bond number which reads Bo = ρgR 2 f /2Γ. The criterion Bo < 1 is always satisfied in our experiments. In a more general prospect, two dimensionless numbers govern the possible motion of a yield-stress fluid drop. The Bingham capillary number Bc that compares yield stress and capillary stress and Bg = σ y /ρgh that compares yield stress and gravity stress and that we call Bingham gravity number in the following. Here h stands for the height of the drop which depends both on the drop size R and on the contact angle θ. To simplify its expression, we can identify h to the initial radius of the drop R 0 and compare both numbers Bc and Bg at the beginning of the spreading. In the limit of small (respectively large) drop size, the Bingham capillary number Bc (respectively Bingham gravity number Bg) is smaller than 1 and the motion is induced by capillary stresses (respectively gravity stresses). Moreover, for sufficiently large yield stresses (σ y > √ Γρg ∼ 25Pa), there is a range of drop sizes, Γ/σ y < R 0 < σ y /ρg, for which motion is impossible. In our experiments, the range of drop radii and heights explored is such as the Bond number is smaller than 1 and the Bingham gravity number is larger than 1. This implies that the spreading of the drop, if it exists, is solely due to capillary forces.

Finally a last point interesting to discuss is the possible role of the elasticity on the spreading of YSF drop.

Depending on the solicitation, a YSF can indeed behave as a visco-plastic fluid or an elastic solid. At the beginning of the spreading, the capillary stress is so high compared to the yield stress that the YSF flows like a liquid [START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF]. But as the spreading goes on, the capillary stress decreases until reaching the yield stress that leads to the drop arrest.

At the end of the spreading, most of the drop behaves like a solid under stress. It is then possible to estimate the stored elastic energy in the solid part of the drop and to compare it to the capillary energy (see Appendix B for more details). Quantitative comparison shows that for Bingham-capillary number B c < 1, and for YSFs with σ y /G < 1, elasticity plays a negligible role. The spreading of a YSF drop allows then to explore mainly its visco-plastic behavior.

V. CONCLUSION

We measure and model the spreading of drops made of YSFs deposited on a wetting solid surface, induced by capillary forces. We observe that after a few minutes, the motion stops and the shape of the drop is frozen. If the Bingham capillary number of the experiment is smaller than 1, the drop takes the shape of a spherical cap. In this limit, the spreading motion stops at a well defined radius and contact angle. We show that the final contact angle does not depend only on the thermodynamical parameters such as interfacial energies but also on the dynamic parameters such as the drop radius, the YSF flow properties and the hydrodynamic boundary condition. With a full analysis of the flow profile in the drop, we predict the contact angle of the arrested state, which corresponds very well to experimental observations. Using a quantitative asymptotic analysis, we show that the final contact angle is directly set by the Bingham capillary number and the hydrodynamic boundary conditions. This macroscopic dynamic arrest then depends on the history of the drop spreading. Here several questions arise concerning the appearance of contact angle hysteresis [START_REF] Eral | Contact angle hysteresis: a review of fundamentals and applications[END_REF], the effect of the impact velocity on the spreading [START_REF] Luu | Drop impact of yield-stress fluids[END_REF], or the existence of a precursor film at the molecular level [START_REF] Popescu | Precursor films in wetting phenomena[END_REF].

FIG. 1

 1 FIG. 1. a) Experimental setup. b) Schematic and notation of a spreading drop: θ the contact angle, R the radius, ξ the height of the drop and zy the thickness of the sheared zone inside the drop. c) Pictures of drops of various YSFs after spreading on rough surface: (I) Carbopol gel (ETD at 0.3%) with a yield stress σy = 6 Pa, (II) and (III) Carbopol gel (ETD at 2%) with σy = 35 Pa, (IV) Carbopol gel (U10 at 0.5%) with σy = 85 Pa. Fits of the drop shape assuming a spherical cap of a same volume as the drop, are plotted in red. Scale bar corresponds to 1 mm.

FIG. 2

 2 FIG. 2. a) Final contact angle θ f versus the final radius R f for three fluids defined in Table I: Carbopol gel (ETD at 2%) with a yield stress σy = 35 Pa (blue diamonds), Carbopol gel (U10 at 0.25%) with a yield stress σy = 29 Pa (red squares) and PAA suspension (at 2%) with no yield stress (green circles), on smooth (open symbols) and rough substrates (plain symbols). b) Final contact angle θ f measured with YSFs on rough substrates as a function of the yield stress σy for a given radius, R f = 1.5 ± 0.2 mm.

FIG. 3 .

 3 FIG. 3. (cos θe -cos θ f ) with θe = 10 • as a function of Bc = R f σy/Γ for two YSFs with various concentrations given in TableI: Carbopol gel (ETD) (blue diamonds) and Carbopol gel (U10) (red squares), on smooth (open symbols) and rough substrates (plain symbols). Plain line corresponds to the numerical resolution assuming a no-slip boundary condition, black dotted line and red dashed line correspond to the analytical asymptotic limit, assuming a no-slip boundary condition (Eq. 3) and a partial slip boundary condition (Eq. 4 with Vs = 0.66V ), respectively. A fitting numerical prefactor of 0.50 ± 0.03 is obtained in both cases. Each data point corresponds to an average over experimental points in a range of 0.1 for the abscissa.

FIG. 4 .

 4 FIG. 4. Schematic of the spreading regimes in the (R f -σy) phase diagram for all YSFs. The dashed line corresponds to the frontier R f σy = Γ. The color scale corresponds to the deviation δ of the final drop shape from a spherical cap as defined in the text (Eq. 5). Roman numbers correspond to the position in this phase diagram of the four pictures from the Fig 1c).

Fig. 4 a

 4 Fig. 4 a phase diagram of all experiments performed as a function of the final radius R f and the yield stress σ y . The color-map corresponds to the value of δ. Pictures of drop profiles (noted I-IV in Fig. 1) are also reported on Fig. 4.

TABLE I .

 I Properties of the complex fluids: type (two Carbopols and one suspension of PAA), weight concentration c, yield stress σy, elastic modulus G and surface tension Γ.

	Type	c (% wt.)	σy (Pa)	G (Pa)	Γ (mN/m)
	Carbopol	3	74 ± 4	206 ± 3	55 * a
	ETD2050	2	35 ± 2	97 ± 1	55 ±10
		1	19.0 ± 1	75 ± 1	58 ±5
		0.3	6.2 ± 0.4	34 ± 0.5	63±3
		0.03	≤1	6.2 ± 0.5	63±2
	Carbopol	0.25	29 ± 2	170 ± 1	55 ± 5
	U10	0.2	18 ± 2	140 ± 1	55 *
		0.15	4.9 ± 0.2	53.0 ± 0.5	55 *
	PAA	2	0	0	60 ± 3

ACKNOWLEDGMENTS

G. Martouzet acknowledges ENS Paris-Saclay for funding. We thank Agnès Piednoir for the measurements of surface roughness by AFM and optical profilometer, Gilles Simon for the conception of the experimental set-up and Yedhir Mezache for interesting discussions.

where c is the dimensionless constant defined in Eq. A10 and V s /V is the fraction of slip at the end of the spreading. This asymptotic limit is plotted in Fig. 3 (see red dotted line with c = 0.5). By imposing V s /V = 0.66, this asymptotic limit describes well the experimental data obtained with Carbopol gels (U10) flowing over smooth surfaces (red open squares) for which a partial slip BC is expected. Concerning the slip velocity, note that we consider a partial slip boundary condition, V s /V < 1, which implies that the slip velocity vanishes when V tends to 0, i.e., when the spreading stops.

Appendix A: Modelisation

In this section, we detail the calculations of the spreading dynamics of the drop until the dynamical arrest is reached and we detail the two asymptotic limits: the beginning and the end of the spreading.

General approach and main assumptions

As mentioned previously, the fluid rheology is assumed to be described by the HB law:

where σ y is the yield stress, K the consistency and n the HB exponent. The geometry of the spreading drop is assumed to be a spherical cap so it is axisymmetric as drawn in Fig. 1b). We also assume that the spreading dynamics and the dynamical arrest are set by a balance between the capillary force and the viscous dissipation. The power injected by the driving capillary force writes P cap = 2πRΓ(cos θ e -cos θ)V , and the power dissipated by viscosity is given by

R 0 ξ a rσ γ dzdr where ξ(r) is the height of the drop and a is a cut-off length that we introduce in order to regularize the divergence of the viscous dissipation at the contact line. In the case of YSFs, the cut-off length is set by the microstructure composing the complex fluids. Here, for Carbopol gels, we assume that the cut-off length is given by the size of a polymer microgel, typically of the order of 1 µm. Balancing the two powers results in:

To calculate the viscous power, it is necessary to determine the shear stress and the velocity everywhere inside the drop during its spreading. The stress results from the momentum balance, and in the lubrication approximation, it reads:

where A is the pressure gradient created by the capillary forces, and ξ is the liquid-air interface position as shown in Fig. 1b). The stress is maximum at the wall and decreases down to zero at the air-liquid interface. From the knowledge of the stress (Eq. A3) and using the HB law (Eq. A1), we can obtain the shear rate everywhere inside the drop:

where z y is the vertical position inside the drop where the shear stress is equal to the yield stress. By definition, σ(z y ) = σ y and z y verifies:

Note that the height of the drop ξ and the pressure gradient A depend both on the radial position r so that z y is also a function of r, making the calculation of the dissipative power more difficult than for a Newtonian fluid. Finally, z y delineates two regions inside the drop: a liquid-like region near the wall and a solid-like region above it. Indeed, for z < z y , γ > 0 and the YSF is sheared and for z > z y , γ = 0, the YSF is not sheared and it flows like a plug, as expected for YSFs [START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF].

We can now integrate the shear rate (Eq. A4) with respect to z and obtain the velocity profile inside the drop:

Here we assume a no-slip boundary condition at the solid surface. Such boundary condition is indeed valid for YSF flowing on rough surfaces [START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF]. From the knowledge of the velocity profile, we can also calculate the mean velocity which depends on the radial position r such as:

In particular, near the contact line, we have lim r→R V m (r) = V .

To solve the spreading dynamics and find R(t), it is necessary first to determine the pressure gradient A(r) by using the mass conservation and then to integrate it in the balance of the powers. This has to be done numerically (see Supplemental Material at URL for details on the numerical resolution). Finally, the link between the angle θ(t) and R(t) is given by the conservation of the drop volume, assuming a spherical-cap shape.

Asymptotic limits

Two asymptotic limits can be solved analytically: 1-the beginning of the spreading characterised by high shear rates and dominated by the viscous properties of the YSF and 2-the end of the spreading characterised by small shear rates and dominated by the plasticity of the YSF. High or low shear rates are defined with respect to γ * = (σ y /K) 1/n for which the two terms of the HB law are equal: σ y = K γn .

At the beginning of the spreading, the shear rate is large compared to γ * such as the yield stress σ y can be neglected with respect to the viscous stress K γn in Eq. A1. So the YSF can be considered as a shear-thinning system without yield stress. The entire problem can be solved analytically and the spreading dynamics follows: R(t) ∝ t p with p = n/(3n + 7), known as the Starov's law [START_REF] Starov | Spreading of non-newtonian liquids over solid substrates[END_REF].

At the end of the spreading, the shear rate is small compared to γ * such as the viscous stress K γn can be neglected with respect to the yield stress σ y . The stress inside the drop is then dominated by the yield stress. Consequently, the power associated to the viscous dissipation can be approximated to: Note that v(ξ) = v(z y ) due the existence of the unsheared region above z y , so that

For YSFs, two boundary conditions (BC) can be encountered: no-slip or partial-slip BCs [START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF][START_REF] Meeker | Slip and flow in soft particle pastes[END_REF]. In the case of a no-slip BC, v(a) = 0 and the dissipative power writes P diss 2πσ y R 0 rv(z y (r)) dr.

To go one step further, we assume, as it is done for Newtonian liquids when deriving the Tanner's law [START_REF] Tanner | The spreading of silicone oil drops on horizontal surfaces[END_REF], that the viscous dissipation dominates in the corner of the drop near the contact line, where regularization is required to avoid divergence of the viscous stress. Using Eq. A7 near the contact line, we can show that the velocity inside the drop and the spreading velocity V are linked by:

The integral over the velocity on the left hand side of Eq. A8 can be calculated knowing the velocity profile inside the drop (Eq. A6):

ξ(r) 0 v r (z) dz = (ξ -z y )v(z y ) + zy 0 v(z) dz. In the asymptotic limit, the sheared region vanishes and becomes negligible compared to the height of the drop, i.e., z y << ξ. This implies that lim σ→σy v(z y ) = V near the contact line. The integral over the radial position in the dissipative power, R 0 rv(z y (r)) dr can be then simplified and approximated to a term proportional to V R 2 . The dissipative power then reads:

where b is a dimensionless constant that takes into account the complexity of the calculation of the integral over the geometry of a spherical cap. The balance of the capillary power and the dissipative power at the end of the spreading leads then to:

where c = b/(2π) is also a dimensionless constant. This asymptotic limit predicts that the distance to the equilibrium (cos θ e -cos θ f ) increases linearly with the Bingham capillary number Bc = (σ y R f )/Γ. This prediction is in very good agreement with what is found experimentally (see Fig. 3). The fit with Eq. A10 of the data, obtained with YSFs spreading on rough surfaces for which no-slip BC is expected leads to c = 0.5 (black dotted line in Fig. 3).

In the case of a partial-slip BC, the velocity v(a) is non zero. By definition of the cut-off length a, it corresponds to the velocity of the first layer of polymer microgel and corresponds to the slip velocity V s which is measured experimentally [START_REF] Péméja | Wall slip regimes in jammed suspensions of soft microgels[END_REF][START_REF] Meeker | Slip and flow in soft particle pastes[END_REF]. Measurements of slip velocity V s showed that it depends essentially on the stress at the wall σ wall , V s (σ wall ) [START_REF] Bonn | Yield stress materials in soft condensed matter[END_REF]. The dissipative power is then:

According to Eq. A3, stress at the wall has the following expression σ wall = A(r)ζ(r) and a priori depends on the radial position r. However, at the end of the spreading, the wall stress tends to the yield stress value and becomes almost independent of r, as shown in Fig. 5. In the asymptotic limit, we can then consider that: σ wall σ y simplifying the calculation of the dissipative power (Eq. A11).

FIG. 5. Stress at the wall σ(z = 0)/σy as a function of the radial position for three velocities: 10 -1 (blue), 10 -3 (orange) and 10 -5 (green). The velocity here is dimensionless, divided by R0(K/σy) 1/n . As the velocity decreases, the wall stress tends to the yield stress. The data are calculated numerically (see Supplemental Material at URL for more details on the numerical resolution).

Finally, as previously done, we assume that the whole dissipation inside the moving drop is dominated by the dissipation inside the corner near the contact line. Using this assumption, we can show that lim σ→σy (v(z y ) -V s (σ wall )) = V -V s . By definition of the dissipative power, V s < V which corresponds to a partial slip boundary condition.

Following the same derivation as in the case of a no-slip BC, we find that:

The balance of powers then leads to:

Appendix B: Role of the elasticity?

We discuss here the possible role of the elasticity (G ) on the spreading of YSF drops. Therefore, we first estimate the elastic energy stored in the arrested drop and compare it to the capillary energy. Second, we show measurements of final contact angles performed with Carbopol gels of same yield stress σ y and different elastic modulus G .

Stored elastic energy

Following the derivation of the velocity profile inside the drop, we show that the drop is composed by a sheared region close to the wall and a non-sheared region for z > z y (see Eq. A6). This non-sheared region behaves like a solid under stress (σ(z y ) = σ y ), which stores elastic energy. At the end of the spreading, the non-sheared region invades the drop and the stored elastic energy is maximum. Assuming a linear elastic response, if we note the local deformation, the elastic energy density is σ = G . Balancing it with the local stress (to replace the value of = σ/G ) and after integration, we can estimate the maximum stored elastic energy:

where V is the volume of the drop. This energy has to be compared to the capillary energy that is bounded by:

Assuming a spherical cap shape, we can show that:

where f depends on the contact angle:

and ranges from 0.2 to 1.2 for 10 • < θ < 70 • .

This estimation shows that elasticity plays a even smaller role as B c and σ y /G are small compared to one.

For typical data as the one shown in Fig. 1c) (picture II) with Carbopol EDT at 2% (σ y = 35 Pa, G = 97 Pa ), a final radius of R f = 1 mm, a final contact angle θ f = 45 • , we find that E elastic /E capillary ∼ 1% which suggests that the stored elastic energy is negligible in the spreading.

Influence of the shear elastic modulus ?

In order to investigate the role of the elastic modulus on the spreading of YSF, we also performed experiments with YSF of similar yield stress values and different elastic modulii. More precisely, we used two carbopol gels (ETD 2050 at 1% and U10 at 0.2%) whose yield stresses are σ y = 18 and 19 Pa respectively and whose elastic moduli are G = 75 and 140 Pa. Although the elastic moduli differ by a factor 2, the final contact angles are identical as shown on Fig. 6 (see blue and red data). These results confirm that the elasticity (G ) has a negligible influence compared to the plasticity (σ y ) on the spreading of YSF drop.

[1] P.-G. de Gennes, F. Brochard-Wyart, and D. Quéré, Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves (Springer, 2004).


Ioan Doré Landau
email: ioan-dore.landau@gipsa-lab.grenoble-inp.fr

Bernard Vau
email: bernard.vau@ixblue.com

Recursive Algorithms for Identification of Dual Youla-Kucera Parametrized Plant Models in Closed Loop Operation

Keywords: System identification, Youla-Kucera parametrization, Recursive algorithms, Identification in closed loop

The growing interest in using dual Youla Kucera plant parametrization for modeling plant uncertainties raises the need for recursive identification algorithms dedicated to the identification of these structures in closed loop operation in view of developing appropriate iterative tuning and adaptive control strategies. The paper presents recursive algorithms for identification in closed loop operation of dual Youla-Kucera parametrized plant models. These algorithms assure global asymptotic stability in the deterministic environment and allow to obtain unbiased parameter estimation in the presence of measurement noise when the plant model is in the model set. The paper also re-visit the Hansen scheme which allows to associate open loop type recursive identification algorithms for the identification of these structures in closed loop operation. When the plant model is not in the model set, comparison of the various algorithms is done in terms of the bias distribution . Further comparisons and performance evaluation is provided by simulations on some relevant examples.

Introduction

There is a growing interest for modeling plant uncertainties using the dual Youla-Kucera plant parametrization [START_REF] Anderson | From Youla-Kučera to identification, adaptive and nonlinear control[END_REF], [START_REF] Mahtout | Advances in youla-kucera parametrization: A review[END_REF] when operating in closed loop. Dual Youla-Kucera parametrization appears to be an efficient way for evaluating robustness of adaptive regulation schemes ( [START_REF] Vau | Adaptive rejection of narrow-band disturbances in the presence of plant uncertainties -a dual youla-kucera approach[END_REF]). Certainly this representation can be used as well for evaluating the robustness of adaptive feedforward disturbance compensation schemes in the presence of plant uncertainties. Application in the field of active noise control have been reported ( [START_REF] Vau | Adaptive rejection of narrow-band disturbances in the presence of plant uncertainties -a dual youla-kucera approach[END_REF]).

Moving towards handling large plant model uncertainties in adaptive regulation and adaptive feedforward disturbance compensation1 will require identification in closed loop of dual Youla Kucera parametrized plant models in order to implement iterative tuning and adaptive control strategies. As far as the authors' knowledge goes, the standard way of recursively identifying a dual parametrized Youla-Kucera plant model in closed loop operation is to use the so called "Hansen scheme" [START_REF] Anderson | From Youla-Kučera to identification, adaptive and nonlinear control[END_REF], [START_REF] Hansen | Closed loop identification via the fractional representation: experiment design[END_REF] which transforms this operation in an open loop type identification scheme. This approach however does not take advantage of using dedicated methods for recursive identification in closed loop operation where the aim is to find a plant model estimate which allows to obtain the best closed loop predictor [START_REF] Landau | An output error recursive algorithm for unbiased identification in closed loop[END_REF], [START_REF] Landau | Recursive algorithms for identification in closed loop. a unified approach and evaluation[END_REF]. The objective is also to get a better approximation in the frequency region close to the Nyquist point when the plant model is not in the model set. There are number of off-line identification procedure for identification of dual Youla Kucera parametrized models in closed loop( [START_REF] Van Den Hopf | Multivariable closed-loop identification: From indirect identification to dual-youla parametrization[END_REF]).

The objective of this paper is to introduce novel recursive identification algorithms dedicated to the identification in closed loop operation of dual Youla Kucera parametrized models which try to minimize the closed loop output prediction error. The algorithms will be developed under the assumptions that the controller is constant and known and that the plant model is in the model set.

The algorithms will assure: global asymptotic stability for any initial parameter estimate and closed loop error prediction, asymptotically optimal prediction for different types of measurement noise, unbiased parameter estimation in the presence of measurement noise under appropriate richness conditions on the excitation signal Then the behavior of the algorithms in the case when the plant model is not in the model set will be analyzed by evaluating the bias in the frequency domain.

The "Hansen scheme" will be re-visited with the objective of enhancing the qualities of the estimated models from the perspective of the bias distribution in the frequency domain when the plant model is not in the model set. A comparison of the various algorithms in terms of bias distribution will be provided.

Simulations on relevant examples will further illustrate the performance of the various algorithms.

The paper is organized as follows: Basic equation will be presented in Section 2.

The algorithms will be developed and analyzed in Sections 3 through 6 under the hypothesis that the plant model is in the model set. The bias distribution in the frequency domain when the plant is not in the model set will be discussed in Section 7. Section 8 will present simulation results .

Basic Equations

The closed loop output error configuration for the identification of a dual Youla-Kucera parametrized plant models is shown in Figure 1 The initial available plant model (nominal model) is described by the transfer operator:

G o (q -1 ) = q -do B o (q -1 ) A o (q -1 ) = B o (q -1 ) A o (q -1 ) (1)
A o = 1 + a 1 q -1 + • • • + a n A q -n A = 1 + q -1 A * 0 ; B o = b 1 q -1 + • • • + b n B q -n B = q -1 B * 0 ; B * o = b 1 + • • • + b n B q -n B +1 ;
where A 0 (q -1 ), B 0 (q -1 ), B * 0 (q -1 ) are polynomials in the delay operator q -1 and n A0 , n B0 and n B -1 represent their orders2 . The feedback fixed polynomial controller K is given by:

K(q -1 ) = R(q -1 ) S(q -1 ) (2)
with:

S = 1 + s 1 q -1 + • • • + s n S q -n S = 1 + q -1 S * ; R = r 0 + • • • + r n R q -n R ;
It will be assumed that S is an asymptotically stable polynomial

The unknown part of the plant model is described by the input output block Π:

Π(q -1 ) = ∆(q -1 ) Γ(q -1 ) (3)
where:

Γ = 1 + γ 1 q -1 + • • • + γ nΓ q -nΓ = 1 + q -1 Γ * ; ∆ = δ 1 q -1 + • • • + δ n∆ q -n∆ = q -1 ∆ * ; (4)
The unknown plant model is described by:

G(q -1 ) = B(q -1 ) A(q -1 ) (5)
and can be expressed as3

G(q -1 ) = ΓB 0 + ∆S ΓA 0 -∆R (6)
The closed loop output predictor can be described similarly except that one uses an estimation of Π:

Π(q -1 ) = ∆(q -1 ) Γ(q -1 ) (7)
where:

Γ = 1 + γ1 q -1 + • • • + γnΓ q -nΓ = 1 + q -1 Γ * ; ∆ = δ1 q -1 + • • • + δn∆ q -n∆ = q -1 ∆ * ; (8)
The input and output of Π will be denoted α(t) and β(t) respectively and the input and output of Π will be denoted α(t) and β(t) respectively.

In the time domain the following input/output relationships will be used:

y(t) = ΓB 0 + ∆S ΓA 0 -∆R u(t) (9)
u(t) = r(t) - R S y(t) (10)
where r(t) is the external excitation.

ŷ(t) = ΓB 0 + ∆S ΓA 0 -∆R û(t) (11)
û(t) = r(t) - R S ŷ(t) (12)
Note that:

α(t) = Su(t) + Ry(t) (13)
β(t) = A 0 y(t) -B 0 u(t) = ∆ Γ α(t) (14) α(t) = S û(t) + Rŷ(t) (15)
β(t) = A 0 ŷ(t) -B 0 û(t) = ∆ Γ α (t) (16)
The closed loop output error is defined as:

ε CL = y(t) -ŷ(t) (17)
and the filtered closed loop output error is defined as:

ν CL = β(t) -β(t) (18)
Using Eqs. ( 10), ( 12),( 14) and ( 16) one effectively gets:

ν CL (t) = P 0 S ε CL (t) (19)
where:

P 0 = A 0 S + B 0 R (20)
is a Hurwitz polynomial. Define:

θ T = [γ 1 , γ 2 . . . δ 1 , δ 2 , . . . ] ; (21) θT = [γ1, γ2 . . . δ1, δ2...], (22)
ϕ T (t) = [-β(t), -β(t -1), . . . α(t), α(t -1), . . . ] (23)
For developing the identification algorithm it will be assumed that the plant model is in the model set. The algorithms will be analyzed first in this context 50 and then in the case when the plant model is not in the model set.

Dual-CLOE Algorithms

In a deterministic context (w = 0 in Fig. 1), for the case of a constant estimated parameter vector θ using the predictor given in Eq [START_REF] Vau | A pseudo-linear regression algorithm in discrete-time[END_REF] one has the following result:

Lemma 1 For a constant estimated vector θ in a deterministic environment (w(t) = 0) the closed loop output error is given by:

ε CL (t + 1) = S ΓP 0 (θ -θ) T ϕ(t) (24)
and the filtered closed loop output error is given by:

ν CL (t + 1) = 1 Γ (θ -θ) T ϕ(t) (25)
The proof of this lemma is given in appendix A.

In the presence of time varying parameter estimates one can define an a priori and and a posteriori filtered closed loop output error, i.e.

ε 0 CL (t + 1) = S ΓP 0 [θ -θ(t)] T ϕ(t) (26)
ε CL (t + 1) = S ΓP 0 [θ -θ(t + 1)] T ϕ(t) (27)
Similarly one can define an a priori and and a posteriori filtered closed loop output error, i.e.

ν 0 CL (t + 1) = 1 Γ [θ -θ(t)] T ϕ(t) (28)
ν CL (t + 1) = 1 Γ [θ -θ(t + 1)] T ϕ(t) (29)
Eqs. ( 27)and (29) have the standard form for using parameter adaptation algorithms (PAA) of the form: [6, 8]

Θ(t + 1) = Θ(t) + F (t)φ(t)ν(t + 1) (30)
ν(t + 1) = ν o (t + 1) 1 + φ T F (t)φ(t) (31)
F (t + 1) -1 = λ 1 F (t) -1 + λ 2 φ(t)φ T (t) (32) 0 < λ 1 < 1 0 ≤ λ 2 < 2, F 0 > 0 4 (33)
In the case of using the closed loop output error one has:

Dual -CLOE algorithm:

Θ(t) = θ(t), φ(t) = ϕ(t) and ν(t) = ε CL (t).
Using the Theorem 4.1 [START_REF] Landau | Recursive algorithms for identification in closed loop. a unified approach and evaluation[END_REF] or Theorem 3.3 [START_REF] Landau | Adaptive control, 2nd Edition[END_REF] the SPR condition by filtering the adaptation error through Γ 0 or filtering the observation vector ϕ(t) through 1/Γ 0 (see next section).

The Dual -FCLOE algorithms

Consider now Eq. ( 24), for the case with time varying estimated parameters.

Neglecting the non-commutativity of the time-varying operators (but an exact algorithms can be derived -see [START_REF] Landau | Adaptive control, 2nd Edition[END_REF]) one can consider to filter the observation vector as follows:

ϕ f (t) = S Γ 0 P 0 ϕ(t) (38)
where Γ 0 is an a priori estimation of Γ (Γ 0 =1 in the absence of a priori information or in the case of FIR models). In this case one gets:

ε CL (t + 1) = Γ 0 Γ [θ -θ(t + 1)] T ϕ f (t) (39)
This equation has the standard form for using the PAA given in Eqs. ( 30) through (33) with:

Dual -F CLOE algorithm: Θ(t) = θ(t), φ(t) = ϕ f (t) and ν(t) = ε CL (t).
The stability condition in this case will be that the transfer function:

Γ 0 (z -1 ) Γ(z -1 ) - λ 2 ; max λ 2 (t) ≤ λ 2 ≤ 2 (40)
be strictly positive real

For the case of using the filtered closed loop output error one can define the regressor vector as:

ϕ f f (t) = 1 Γ 0 ϕ(t) (41)
Eq. ( 25) becomes:

ν CL (t + 1) = Γ 0 Γ (θ -θ(t + 1)) T ϕ f f (t) (42)
In this case, the algorithm will take the form: estimation of Γ(t) but including a stability test (Dual-AFCLOE algorithm). .

The Noisy Case

In the presence of measurement noise w(t) , the output of the plant will be given by:

y((t + 1) = -Ay(t) + Bu(t) + Aw(t + 1) (43)
One assumes that the measurement noise w(t) is bounded and independent with respect to the external excitation r(t) and that one uses a decreasing adaptation gain (λ 1 (t) ≡ 1, λ 2 (t) > 0). One can use for analysis, the ODE approach [9].

Following a similar path as for the proof of Lemma 1 it can be shown that the closed loop output prediction error and the filtered closed loop output error for a fixed value of the estimated parameter vector will be given by:

ε CL (t + 1) = S ΓP 0 (θ -θ) T ϕ(t) + AS ΓP 0 w(t + 1) (44) ν CL (t + 1) = 1 Γ (θ -θ) T ϕ(t) + A Γ w(t + 1) (45)
For the Dual-FCLOE algorithm the expression of the adaptation error for a constant estimated paramter vector will take the form:

ε CL (t + 1) = Γ 0 Γ (θ -θ) T ϕ f (t) + AS ΓP 0 w(t + 1) (46)
ν CL (t + 1) = Γ 0 Γ (θ -θ) T ϕ f f (t) + A Γ w(t + 1) (47)
One observes that ϕ(t), ϕ f (t) and ϕ f f (t) are independent with respect to w(t+1)

(by the way that they are generated) and one can straightforwardly use Theorem 4.1 [START_REF] Landau | Adaptive control, 2nd Edition[END_REF] leading to the following conclusions:

80
Assume that the stationary process φ(t, θ) and ν(t + 1, θ) can be defined for θ(t) ≡ θ.

Assume that θ(t) generated by the algorithm belongs infinitely often to the domain (D s ) for which the stationary processes φ(t, θ) and ν(t + 1, θ) 85 can be defined.

Assume that either w(t + 1, θ) is a sequence of independent equally distributed normal random variables (0, s), or E{φ(t, θ), w(t + 1, θ)} = 0

Then, if the strictly positive real condition of Eq. ( 35) for Dual-CLOE, of Eq. (37) for FDual-CLOE, of Eq. (40) Dual-FCLOE and FDual-FCLOE are satisfied, one has:

P rob{ lim t→∞ θ(t) ∈ D c } = 1 (48)
where the convergence domain D c is defined as :

D c : θ : φ T (t, θ)[θ -θ] = 0 (49)
Under richness conditions the domain D c reduces to θ = θ and one gets unbiased parameter estimates i.e.

P rob{ lim

t→∞ θ(t) = θ} = 1 (50)

Dual -XCLOE Algorithm

In order to remove the SPR condition (37) depending upon the unknown polynomial Γ(q -1 ) one can consider to use an extended estimation model. However, the independence between the measurement noise and the regressor vector will be lost in the noisy case even for constant estimated parameters. So in order to remove the SPR condition (37) and to assure unbiased estimates in the presence of noise, we will consider an ARMAX model representation for the plant + noise:

y((t + 1) = -Ay(t) + Bu(t) + Ce(t + 1) (51)
where e(t) is a gaussian discrete time white noise (0,σ) and

C = 1 + q -1 C * is
an asymptotically stable polynomial. The objective will be first to construct an optimal predictor for the case of known parameters which will assure that the prediction error will become asymptotically white and then an adaptation algorithm will be developed for the case of unknown parameters. One has the following result:

Lemma 3 Under the hypotheses that A 0 , B 0 , R, S, ∆, Γ and C are perfectly known, the asymptotically optimal one step ahead closed loop output predictor is given by:

Γ β(t + 1) = ∆ * α(t) + Ĥ * ε CL S (52) ŷ(t + 1) = -A * 0 ŷ(t) + β(t + 1) + B * 0 û(t) (53)
where: Note that Eq. ( 52) can be rewritten as:

H * = C * S -(ΓA 0 ) * S -Γ * B * 0 R : H = 1 + q -1 H * Ĥ * = ĥ1 + ĥ2 q -1 + ... + ĥn H q -n H +1 ; Ĥ = 1 + q -1 Ĥ * n H = max(n C + n S , n Γ + n A0 + n S , n Γ + n B0 + n R ); (54)
β(t + 1) = θT e ϕ e (t) (56)
where:

θT e = [ θT , ĥ1 , ĥ2 , ..., ĥnH ] (57) ϕ T e (t) = [ϕ T , ε CLf (t), ...ε CLf (t -n H + 1)] (58) ε CLf (t) = 1 S ε CL (t) (59)
To identify the vector of parameter, one replaces θe in the predictor by a timevarying vector θe (t) which will be estimated.

In a deterministic context e(t + 1) = 0 and C = 1. From Eqs. (99) through (102) it results that the a posteriori output of the Π predictor will be given by:

β(t + 1) = θT e (t + 1)ϕ e (t) (60)
and the a posteriori closed loop prediction error will be given by:

ε CL (t + 1) = [θ e -θe (t + 1)] T ϕ e (t) (61)
This equation has the standard form leading to the use of the parameter adaptation algorithm given in Eqs. (30) through (33) with;

Dual-XCLOE algorithm:

Θ(t) = θe , φ(t) = ϕ e (t), ν(t) = ε CL (t).
Lemma 3bis Under the hypotheses that A 0 , B 0 , R, S, ∆, Γ and C are perfectly known, the asymptotically optimal one step ahead filtered closed loop output predictor is given by:

Γ β(t + 1) = ∆ * α(t) + Ĥ * ν CL (62)
where:

H * = C * -Γ * ; Ĥ * = C * -Γ * ; Ĥ = 1 + q -1 Ĥ * = ĥ1 + ĥ2 q -1 + ... + ĥn H q -n H +1 ; n H = max(n C , n Γ ); (63)
leading to:

lim t→∞ [ν CL (t + 1) -e(t + 1)] = 0 ( 64
)
for θ = θ and Ĥ = H. The proof of this lemma is given in Appendix B. Note that Eq. ( 62) can be rewritten as:

β(t + 1) = θT ef ϕ ef (t) (65)
where:

θT ef = [ θT , ĥ1 , ĥ2 , ..., ĥn H ] (66)
ϕ T ef (t) = [ϕ T , ν CL (t), ...ν CL (t -n H + 1)] (67)
To identify the vector of parameters, one replaces θef in the predictor by a time-varying vector θef (t) which will be estimated.

In a deterministic context e(t + 1) = 0 and C = 1, the a posteriori output of the predictor will be given by:

β(t + 1) = θT e (t + 1)ϕ ef (t) (68)
and the a posteriori prediction error will be given by:

ν CL (t + 1) = [θ ef -θ ef (t + 1)] T ϕ ef (t) (69)
This equation has the standard form leading to the use of the parameter adap-90 tation algorithm given in Eqs. (30) through (33) with:

FDual-XCLOE algorithm:

Θ(t) = θef , φ(t) = ϕ ef (t), ν(t) = ν CL (t).
Remark: The algorithm using the adjustable predictor for the filtered closed loop output error will need less parameters to estimate than the algorithm us-95 ing the adjustable predictor for the closed loop output error.

Stability analysis

Using the Theorem 4.1 [START_REF] Landau | Recursive algorithms for identification in closed loop. a unified approach and evaluation[END_REF] or Theorem 3.3 [START_REF] Landau | Adaptive control, 2nd Edition[END_REF] together with the boundedness of ν CL and ϕ e (t) for any initial conditions.

The noisy case

In the noisy case, using Eqs. (99) and (100)from Appendix B, the prediction error equation takes the form:

ε CL (t + 1) = [θ e -θe (t + 1)] T ϕ e (t) + Ce(t + 1) -C * ε CL (t) (72)
which leads to:

ε CL (t + 1) = 1 C [θ e -θe (t + 1)] T ϕ e (t) + e(t + 1) (73)
One can use for analysis the ODE method [9]. Since e(t+1) is a white noise sequence, ϕ e ( θe , t) will be independent with respect to e(t+1) and applying Theorem 4.1 from [START_REF] Landau | Adaptive control, 2nd Edition[END_REF] one concludes that asymptotic unbiased parameter estimates are obtained under richness conditions provided that

1 C - λ 2 2 , λ 2 = maxλ 2 (t) (74)
is a strictly positive real transfer function. The same convergence result is ob-100 tained when using the predictor of Eq. (62) (filtered closed loop output predictor) with adjustable parameters and the corresponding parameter and regressor vector.

A summary of the algorithms presented above is given in Table 1.

Modified Hansen approach for recursive identification of Youla

l(t) = -B 0 u(t) + A 0 y(t) (75b)
and on the other hand

α(t) = Su(t) + R(y(t) -w(t)) (76a)
β(t) = -B 0 u(t) + A 0 (y(t) -w(t)) (76b)
But since

β(t) = ∆ Γ α(t) (77)
and combining with ( 5) and ( 6), one obtains:

l(t) = ∆ Γ p(t) + A Γ w(t) (78)
If one uses F-OLOE algorithm [START_REF] Landau | Adaptive control, 2nd Edition[END_REF] in order to identify the Dual-YK filter with this scheme, the parameters vectors remain the one given by ( 21) and ( 22), the observation vector is ϕ

T (t) = [-l(t), -l(t -1) • • • p(t), p(t -1) • • • ], the
predicted output is l(t + 1) = θT (t)ϕ(t), the regressor vector is ϕ(t) filtered by 1 Γ0 , and the adaptation error is given by ν(t + 1) = ε(t + 1) = l(t + 1)l(t + 1) (MH-FOLOE algorithm). Another choice consists in using an ARMAX predictor model, and one can use the recursive extended least-squares (RELSsee [START_REF] Landau | Adaptive control, 2nd Edition[END_REF]), where the regressor and the observation vector are in this case: ϕ T e (t) = [ϕ T (t), ε(t), ε(t -1) • • • ], and the parameter vector is θT e = [ θT , ĉ1 , ĉ2 • • • ]. This structure can be compared to the well-known Hansen scheme where an open-loop identification is performed with input-output sequences {x(t)} and {z(t)} defined as

x(t) = S 2 P 0 r(t) (79a) z(t) = y(t) - B 0 S P 0 r(t) (79b)
One can verify that

z(t) = ∆ Γ x(t) + AS ΓP 0 w(t) (80)
Consequently the Hansen structure induces an implicit frequency weighting by

S 2 (e iω )
P0(e iω ) with respect to the excitation sequence {r(t)} which in general has not a specific interest, whereas the scheme proposed in this section is a truly open-loop procedure (without any weighting) for the identification of the the Dual Youla-Kucera model in closed loop operation.

Alg. Dual-CLOE Dual-FCLOE Dual-XCLOE FDual-CLOE FDual-FCLOE FDual-XCLOE Plant+Noise y = B A + w y = B A + w y = B A + C A e y = B A + w y = B A + w y = B A + C A e Adj. Par. Vect. Θ θT = [γ 1 . δ1 .] θT == [γ 1 . δ1 .] θT e = [γ 1 . δ1 . ĥ1 .] θT = [γ 1 . δ1 .] θT = [γ 1 . δ1 .] θT e = [γ 1 . δ1 . ĥ1 .] Obs. vect. ϕ [-β(i), α(i)] [-β(i), α(i)] [-β(i), α(i), ε f (i)] [-β(i), α(i)] [-β(i), α(i)] [-β(i), α(i), ν f (i)]
Pred.out. ŷ0 (t+ 1)

-A * 0 y + B * 0 û + β0 -A * 0 y + B * 0 û + β0 -A * 0 y + B * 0 û + β0 β = θT (t)ϕ(t) β = θT (t)ϕ(t) β = θT (t)ϕ(t) Adapt err. y -ŷ y -ŷ y -ŷ β -β β -β β -β Regr.vect. φ ϕ S Γ 0 P 0 ϕ ϕ ϕ S Γ 0 P 0 ϕ ϕ Stab.cond. S ΓP 0 - λ 2 2 Γ 0 Γ - λ 2 2 None 1 Γ - λ 2 2 Γ 0 Γ - λ 2 2 None (deterministic) SPR SPR SPR SPR Conv. cond.) S ΓP 0 - λ 2 2 Γ 0 Γ - λ 2 2 1 C - λ 2 2 1 Γ - λ 2 2 Γ 0 Γ - λ 2 2 1 C - λ 2

Asymptotic bias distribution in the frequency domain

It is important to assess the bias distribution in the frequency domain of the estimated model when the plant model is not in the model set. Contrary to prediction error methods which aim at minimizing the prediction error variance asymptotically, it has been shown in [START_REF] Vau | Some remarks on the bias distribution analysis of discrete-time identification algorithms based on pseudo-linear regressions[END_REF], that in general this is not the case for the algorithms belonging to the Pseudo-Linear Regression (PLR) class, like Dual-CLOE. The non-measurable signal whose variance is asymptotically minimized by algorithms of the PLR class, is called the equivalent prediction error denoted as ε E , and from the expression of this latter signal it is possible to infer a limit model in the frequency domain, useful to analyze the bias distribution of the algorithm. The proof is given in Appendix D (section 14). where Ψ p (ω) is the PSD of the excitation p, when the identification of the Dual-YK filter is performed with F-OLOE [START_REF] Landau | Adaptive control, 2nd Edition[END_REF] (where the filter of the observation vector is 1/Γ 0 ), or with H-OLOE [START_REF] Vau | A pseudo-linear regression algorithm in discrete-time[END_REF] (where the basis poles are the roots of Γ 0 ). The proof is given in Appendix F (section 16).

By comparison, the limit model of the Hansen scheme, when identified with F-OLOE where the observation vector is the regressor filtered through 1 Γ0 , or with H-OLOE where the basis poles are the roots of Γ 0 (it can be inferred from the one of PEM [9]), is given by: θ * = Argmin 2 Ψ r (ω)dω (86)

Comparison of the various algorithms in terms of bias

The Dual-FCLOE algorithm seems to offer the best approximation properties since the model errors in the frequency domain are weighted by the true and the estimated sensitivity function. This means that the best approximation is obtained in the the frequency zone close the Nyquist point (where the modulus of the output sensitivity function has its maximum). The other algorithms from the CLOE family induce a larger penalty in high frequencies.

Simulations

One considers the following configuration:

B 0 (q -1 ) = q -1 + 0.5q -2 A 0 (q -1 ) = 1 -1.5q -1 + 0.7q -2 R(q -1 ) = 0.2419 -0.1491q -1 -0.2144q -2 + 0.1766q -3 S(q -1 ) = 1 -0.5043q -1 -0.3696q -2 -0.1262q -3 ∆(q -1 ) = 2.8q -1 -0.7q -2 Γ(q -1 ) = 1 -1.1q -1 + 0.3q -2
The simulations are performed under the following conditions:

The excitation is a PRBS (2048 samples, decimation 1). The sequence {w(t)} is a centered gaussian white noise, and in closed-loop the noise/noise ratio is about 32% (-10dB) (output variance/noise variance) for all simulations 5 . In this context we are far from the asymptotic behaviour of the algorithms but closer to real situations. Fig. 3 displays the Bode Diagrams of the nominal system (B 0 , A 0 ), and the current one (note that this latter is unstable in open-loop, but is stabilized by the controller), and Fig. 4 compares the output sensitivity functions of the loops including these two systems. Note that the maximum of the sensitivity function including the true system is roughly equal to 10 dB, which implies a modulus margin lower than 0.32. The first set of simulations compares 50 realizations in the case "plant model is in the model set" for the following algorithms:

a) The Hansen scheme where the identification routine is F-OLOE (or equivalently H-OLOE see [START_REF] Vau | A pseudo-linear regression algorithm in discrete-time[END_REF]), where the zeros of the filter (or the basis poles of H-OLOE) are equal to 0.5 (Fig. 5),

b) The modified Hansen scheme, where the parameter estimation is performed with F-OLOE or H-OLOE, and where the zeros of the filter (or the basis poles of H-OLOE) are equal to 0.5 (Fig. 6), c) Dual-FCLOE, where the roots of Γ 0 are set to 0.5 (Fig. 7), d) FDual-FCLOE, where the roots of Γ 0 are set to 0.5 (Fig 8).

These figures reveal that in this highly noisy context the Hansen scheme provides the lowest variance in the low and high frequency ranges and the modified

Hansen Scheme is the one yielding the larger variance in the low and high frequency range. However these frequencies domains are not essential for control.

In the critical frequency range (around the maximum of the sensitivity function and of the system resonance) all algorithms show a very little difference One can notice that the misfit is strongly penalized for the frequencies around the maximum of the sensitivity function. As expected from expression (83) , F Dual -CLOE penalizes more significantly the fit in high frequency (Fig. 11), whereas the modified Hansen scheme, in accordance with Eq. (85) spreads the

Validation of the identified models

The techniques presented in [START_REF] Landau | Recursive algorithms for identification in closed loop. a unified approach and evaluation[END_REF] and [START_REF] Landau | Adaptive control, 2nd Edition[END_REF] can be used also in this context.

Conclusions

The paper has presented a set of recursive algorithms for identification in closed loop operation of plant models represented in dual Youla Kucera form.

The algorithms have been analyzed, compared and evaluated by simulations. In this context it appeared that there are three basic options for the choice of the prediction error:1) the closed loop output error 2) the closed loop error between the outputs of the true and estimated dual YK block which corresponds to a filtered closed loop output error) and 3) The Hansen scheme (and its modification) combined with open loop recursive identification algorithms. Some of these algorithms fully exploit the fact that one operates in closed loop allowing to get identified models which provide a better approximation in the crucial frequency regions close to the Nyquist point since the bias is heavily weighted by the sensitivity functions. Some of the algorithms allow to get models without any positive real condition related to a priori estimates.

12. Appendix B: Proof of Lemmas 3 and 3bis

From Eq. (51) using the expression of A and B given in Eq. ( 6) one gets:

Γ[A 0 y(t + 1 -B 0 u(t + 1] = ∆[Su(t + 1) + Ry(t + 1] + Ce(t + 1) (95)
which can be alternatively written as:

y(t + 1) = -(ΓA 0 ) * y(t) + ΓB 0 u(t + 1) +∆ * [Su(t) + Ry(t)] + Ce(t + 1) (96)
where:

ΓA 0 = 1 + q -1 (ΓA 0 ) * (97)
Adding and subtracting -(ΓA 0 ) * ŷ(t) + ΓB 0 û(t + 1) + ∆ * [S û(t) + Rŷ(t)] one gets:

y(t + 1) = -(ΓA 0 ) * ε CL (t) -ΓB * 0 R S ε CL (t) -(ΓA 0 ) * ŷ(t) + ΓB 0 û(t + 1) + ∆ * [S û(t) + Rŷ(t)] + Ce(t + 1) +∆ * [Rε CL (t) -S R S ε CL (t)] (98)
where the last parenthesis is null. Adding and subtracting the term C * ε CL (t),

Eq. (98) will take the form:

y(t + 1) = [C * S -(ΓA 0 ) * S -Γ * B * 0 R] ε CL (t) S -(ΓA 0 ) * ŷ(t) + ΓB * 0 û(t) + ∆ * [S û(t) + Rŷ(t)] -C * [ε CL (t) -e(t)] + e(t + 1) (99)
From Eq. (52) one gets:

ŷ(t + 1) = -( ΓA 0 ) * ŷ(t) + ΓB * 0 û(t) + ∆ * [S û(t) + Rŷ(t)] + Ĥ * ε CL S (100)
Subtracting now Eq. (100) from Eq.(99) for θ = θ and taking into account that:

(ΓA 0 ) * -( ΓA 0 ) * = (Γ * -Γ)A 0 (101) (Γ -Γ) = q -1 (Γ * -Γ * ) (102)
yields:

ε CL (t + 1) = -C * (ε CL (t) -e(t) + e(t + 1) (103)
which implies that (55) holds since C(z -1 ) has all his roots inside the unit circle. Proof of lemma 3bis

Eq.( 95) can be rewritten using Eqs.( 13) and ( 14):

β(t + 1) = -Γ * β(t) + ∆ * α(t) + Ce(t + 1) = -Γ * ν CL (t) -Γ * β(t) + ∆ * α(t) + Ce(t + 1) (104)
by adding and subtracting the term -Γ * β(t)+∆ * α(t). Eq. ( 62) can be rewritten Case 1: ν(t) = ε CL (t). As in [START_REF] Vau | Some remarks on the bias distribution analysis of discrete-time identification algorithms based on pseudo-linear regressions[END_REF], let us define Q(q -1 , θ) such that Q(q -1 , θ) ∂ε CL ∂ θ = -ϕ(t, θ). Notice that we have β(t + 1) = P0 S ŷ(t + 1, θ) -B 0 r(t), and β(t + 1) = θT ϕ(t). Therefore ∂ ŷ(t+1,∂ θ) Case 2: ν(t) = ν CL (t). The expression of Q(q -1 , θ) such that Q(q -1 , θ) ∂ν CL ∂ θ = -ϕ(t, θ) leads to Q(q -1 , θ) = Γ, since β(t + 1) = θT ϕ(t, θ). One has ν CL (t) = P0 S ε CL (t), and the equivalent prediction error ε E (t) is given by ε E (t) = Q(q -1 , θ)ν CL (t)+ (1 -Q(q -1 , θ)) AS ΓP0 w(t + 1), and one gets ε E (t) = Â G -Ĝ AS ΓP0 r(t) + AS ΓP0 w(t), leading to the same limit model as in case 1.

Appendix D: Proof of Lemma 5

Like in case 1 of Dual -CLOE (where ν(t) = ε CL (t)), one has Q(q -1 , θ) = P0 Γ S . On the other hand, one has ϕ f (t) = S P0Γ0 ϕ(t) = 1 Q f (q -1 ) ϕ(t), where Q f (q -1 ) = P0Γ0 S . Therefore one can apply Theorem 2 in [START_REF] Vau | Some remarks on the bias distribution analysis of discrete-time identification algorithms based on pseudo-linear regressions[END_REF], with the present expressions of Q(q -1 , θ) and Q f (q -1 ), and the equivalent prediction error is ε E (t) = Q(q -1 , θ)ν CL (t) + (1 -Q(q -1 , θ)) AS Γ0P0 w(t + 1), and one gets ε E (t) = ÂS Γ0P0 G -Ĝ AS ΓP0 r(t)+ AS ΓP0 w(t), hence the limit model in the frequency domain.

Appendix E: Proof of Lemma 7

Let us determine Q(q -1 , θ) such that Q(q -1 , θ)

Appendix F: Proof of Lemma 8

One has l(t) = ∆ Γ p(t) + A Γ w(t). Let us define l(t) as the predicted output of F-OLOE (or H-OLOE), one has l(t) = ∆ Γ p(t). Set ν(t) the prediction error of F-OLOE (or H-OLOE) such that ν(t) = l(t) -l(t). On the other hand l(t) = -B 0 u(t) + A 0 y(t) and u(t) = -R S y(t) -1 S p(t). Therefore one obtains ν(t) = P0 S ε(t), where ε(t) = y(t) -ŷ(t). Moreover, one has r(t) = -1 S p(t). Thus

Figure 1 :

 1 Figure 1: Closed loop ouput error identification of dual Youla Kucera parametrized model

  leading to: lim t→∞ [ CL (t + 1) -e(t + 1)] = 0 (55) for θ = θ and Ĥ = H. The proof of this lemma is given in Appendix B.

Figure 2 :

 2 Figure 2: Dual YK identification from a YK structure (without Q-filter)

105

  Kucera structure in closed loop operationAdding a Youla-Kucera observer whose output is {l(t)} (it is the input to a Youla-Kucera control filter), a direct identification of the Dual-Youla Kucera structure is possible by introducing an excitation sequence {p(t)} and collecting the sequence {l(t)} as shown in Fig.2. The identification can be performed with any open-loop type algorithm. One has indeed p(t) = Su(t) + Ry(t) (75a)

7. 1 .Lemma 4 Â(e iω ) 2 G 2 S

 1422 Limit model for the Dual-CLOE algorithm One has the following result 125 The limit model for Dual -CLOE algorithm ( ν(t) = ε CL (t)) or for the F Dual -CLOE algorithm ( ν(t) = ν CL (t)) is given by θ * = Argmin +π -π (e iω ) -Ĝ(e iω ) yw (e iω ) 2 Ψ r (ω)dω (81)where S yw = AS ΓP 0 is the output sensitivity function and Ψ r (ω) is the PSD (Power Spectral Density) of the external excitation signal.The proof is given in Appendix C (section13).

7. 2 .Lemma 5 2 G 2 S

 2522 Limit model for the Dual-FCLOE algorithm 130 The limit model of Dual-FCLOE algorithm is given by θ * = Argmin +π -πÂ(e iω )S(e iω ) Γ 0 (e iω )P 0 (e iω ) (e iω ) -Ĝ(e iω ) yw (e iω ) 2 Ψ r (ω)dω (82)

7. 3 . 6 2 G 2 S

 3622 Limit model for the F Dual-FCLOE algorithm Lemma The limit model of F Dual-FCLOE algorithm is given by θ (e iω ) -Ĝ(e iω ) yw (e iω ) 2 Ψ r (ω)dω (83) The proof is obtained by combining results of Appendices C and D (sections 13 and 14).

Lemma 7 Â(e iω ) 2 G 2 SÂ(e iω ) 2 WÂ(e iω ) 2 Ψ 8 2 G 2 A 2 Ψ

 722228222 The limit model of Dual-XCLOE algorithm is given by θ * = Argmin +π -π (e iω ) -Ĝ(e iω ) yw (e iω ) 2 Ψ r (ω)+ (e iω ) A(e iω ) Γ(e iω ) Â(e iω )Γ(e iω ) -Ĉ(e iω ) e (ω)dω (84) where Ψ e (ω) is the PSD of the noise The proof is given in Appendix E (section 15) 7.5. Limit model for the MH-FOLOE Lemma The limit model of this algorithm is given by θ * = Argmin +π -π Â(e iω ) Γ 0 (e iω ) (e iω ) -Ĝ(e iω ) (e iω ) P 0 (e iω )Γ(e iω ) p (ω)dω (85)

2 G 2 S

 22 iω )S(e iω ) Γ 0 (e iω )P 0 (e iω ) (e iω ) -Ĝ(e iω ) yw (e iω

  )

Figure 3 :

 3 Figure 3: Bode diagrams of the nominal system and the current one

Figure 4 :Figure 5 :

 45 Figure 4: Sensitivity functions for the closed-loop including the nominal system and the current one

Figure 6 :Figure 7 :Figure 8 :j=1 σ 2

 6782 Figure 6: Identification result (50 realizations)for the modified Hansen scheme, estimation performed with F-OLOE (MH-OLOE)

Figure 9 :

 9 Figure 9: Resulting model for a Dual-YK Filter of order one, Hansen scheme

Figure 10 :Figure 11 :

 1011 Figure 10: Resulting model for a Dual-YK Filter of order one, Dual-FCLOE

  + 1) = -Γ * β(t) + ∆ * α(t) + [-Γ * + C * ]ν CL (t) (105)Subtracting Eq. (105) from Eq. (104) one gets for the known parameter case:ν CL (t + 1) = -C * (ν CL (t) -e(t) + e(t + 1)(106)which implies that Eq. (64) holds.13. Appendix C: Proof of Lemma 4225

  t). Therefore Q(q -1 , θ) = P0 S + θT ∂ϕ qε = P0 Γ S . Now, Theorem 1 of [10] can be used with the present expression ofQ(q -1 , θ) = P0 Γ S . Since from [8], 230 p.308, ε CL (t) = ÂS P0 Γ G -Ĝ AS ΓP0 r(t) + AS ΓP0 w(t).One obtains, the equivalent prediction error ε E (t) = Â G -Ĝ AS ΓP0 r(t) + AS ΓP0 w(t), and one has θ * = Argmin +π -π Ψ ε E (ω)dω, where Ψ ε E (ω) is the PSD of ε E , hence the limit model in the frequency domain.

θ

  ∂ε CL ∂ θ = -ϕ(t, θ). One has β(t + 1) = θT e φ e (t), and ŷ(t + 1) = S P0θT e φ e (t). Therefore ∂ ŷ(t+1) , andQ(q -1 , θ) = P0 S + θ ∂ϕe ∂qε = P0 S + ε E (t) = Â (G -Ĝ) AS P0Γ r(t) + W A Γ ÂΓ -Ĉ Â e(t), like in the classical X-CLOE algorithm.

Table 1 :

 1 Summary of closed loop identification algorithms for dual Youla Kucera plant models

	2

Table 2 :

 2 Standard deviation of the estimated parameters

	True	MH-OLOE	Dual-	FDual-	Hansen
	Parameters		FCLOE	FCLOE	Scheme/F-
					OLOE
	γ 1 = -1	0.1939	0.1121	0.1981	0.1137
	γ 2 = 0.3	0.1393	0.0783	0.1443	0.0820
	δ 1 = 2.8	0.1569	0.0880	0.0965	0.0580
	δ 2 = -0.7	0.5938	0.4096	0.6602	0.3894
	190
	The second set of simulations consists in performing identification in the case
	where the true Dual-Youla parameter is not in the model set (here ∆ and	Γ

Table 3 :

 3 Standard indicators of the parameters' dispersion

These schemes are adaptive with respect to variations of dynamic characteristics of the disturbances but the plant model is assumed to be known and almost invariant

The complex variable z -1 will be used for characterizing the system's behavior in the frequency domain and the delay operator q -1 will be used for describing the system's behavior in the time domain.

In order to simplify the writing the argument q -1 has been dropped out in many equations

λ 1 (t) and λ 2 (t) allow to obtain various profiles for the evolution of the adaptation gainF (t)

It corresponds to a high level of the measurement noise compared to the usual ratio used in the literature which is around 10% (-20 dB)

In the case of using the filtered closed loop output error one has: F Dual -CLOE algorithm: Θ(t) = θ(t), φ(t) = ϕ(t) and ν(t) = ν CL (t).

Stability Analysis

FDual -FCLOE algorithm:

and the stability condition will be that the transfer function given in Eq. ( 40) is a strictly positive real (SPR) transfer function.

Remark: One can also consider to filter the regressor vector through the current bias over the all spectrum, as shown in Fig. 12.

Appendix A: Proof of Lemma 1

Proof : From Eq. ( 14) one has

which can be alternatively written as:

Adding and subtracting the term -Γ * β(t) + ∆ * α(t) one gets:

From Eq. ( 16) one gets

Subtracting Eq. (90) from Eq. ( 89) and using Eq. ( 18) one gets:

By using Eqs.( 10),( 12),(13),(15) one concludes that:

and using now Eqs. ( 21) through (23), it results that Eq. ( 91) takes the form:

from which Eq. ( 25) results. Using now Eqs. ( 18),( 14),( 16), ( 10),( 12),(17), one gets

Combining Eqs. ( 93) and (94) one gets Eq. ( 24)

End of the proof.

one has ν(t) = Â Γ G -Ĝ A P0Γ p(t) + A Γ w(t). Now the equivalent prediction error ν E (t) of F-OLOE or H-OLOE is given by ν

And finally ν E (t) = Â γ0 (G -Ĝ) A P0Γ p(t) + A Γ w(t), hence the limit model in the frequency domain, and this ends the proof.

for the efficient identification of stiff systems, in: Proceedings of the 19th IFAC symposium on system identification (Sysid 2021), 2021.