G Reinisch
email: reinisch@lcts.u-bordeaux1.fr

J.-M Leyssale

N Bertrand

G Chollon

F Langlais

G Vignoles

Experimental and theoretical investigation of BCl 3 decomposition in H 2

Keywords: Chemical Vapor Deposition, Boron, Boron Carbide, Ab-initio computation, Chemical Kinetics

A combined experimental and theoretical study of the homogeneous decomposition of BCl 3 in a H 2 carrier gas is presented. A detailed description of the B/Cl/H thermodynamic equilibrium is first obtained from ab initio calculations from which a restricted low energy chemical mechanism is identified to model the decomposition of BCl 3 . Transition state theory is then invoked to obtain reaction rates and the resulting kinetic mechanism is incorporated in a 1D model of a CVD reactor. Comparison of calculated steady state concentrations with in situ FT-IR measurements shows a good agreement at low temperatures, thus validating the kinetic model. The divergence observed at higher temperatures is attributed to boron deposition.

Introduction:

Self-healing Ceramic-Matrix Composites (CMCs) are promising materials for demanding aeronautic applications, i.e. simultaneous exposition to heat, mechanical loads, and oxidizing atmospheres. Their internal protection against oxidation can be based on boron-containing matrix layers, because of their ability to generate a liquid oxide, which acts as a diffusion barrier [1]. Accordingly, the preparation of these materials relies on the precise control of the boron-containing layer thicknesses. Chemical Vapor Deposition (CVD) is the privileged method for their production, since it allows, by a precise tuning of the reaction to diffusion ratio, to achieve deposit uniformity inside a complex fibrous medium. In that case the process is termed CVI (Chemical Vapor Infiltration) and operates usually at low pressures and moderate temperatures.

As boron hydrides like BH 3 are extremely toxic and hazardous, boron halides and in particular BCl 3 are commonly used as boron precursors. The decomposition of BCl 3 in a hydrogen-rich CVD reactor is however not fully characterized. Most experimental papers [2][3][4][5] report only on deposition rates, while little or no information is given on the gas phase decomposition of BCl 3 . Nonetheless, Sezgi et al. [2] indicate that BHCl 2 is the main gaseous product in atmospheric pressure CVD. Several modeling works have been performed in parallel to these experimental characterizations; featuring separately or together: thermochemical equilibrium computations, mass-transfer modeling, and chemical kinetic mechanisms. Two gas phase mechanisms, based on ab-initio calculations, have been proposed. Harris et al. [6] proposed a mechanism based on 31 chemical reactions. The decomposition of BCl 3 starts with the reaction BCl 3 +H→BCl 2 +HCl, which depends itself on a first H 2 dissociation step. No direct validation with experimental data has been given. More recently, in the frame of a B-doped Si CVD modeling study, Lengyel and Jensen have proposed a restricted mechanism [7], which takes only into account the three bimolecular reactions listed later in table 2. None of these reactions were involved in the previous mechanism [6]. No direct validation with respect to gas-phase experimental data has been given; however the boron mass uptake data have been reproduced using fitted sticking coefficients.

Accordingly, there are still some open questions concerning the decomposition of BCl 3 in H 2 , especially at low pressures; this has motivated the present study. Here, the emphasis is set on the gas phase characterization and modeling, from ab-initio methods to reactor simulation.

Experimental:

Decomposition experiments are performed in a low pressure hot wall tubular reactor. A mixture of BCl 3 and H 2 flows through a vertical silica tube heated by a graphite susceptor coupled to a radio-frequency inductive coil. The tube section is 34 mm in diameter and the hot isothermal zone length is 140 mm. The reactor outlet is connected to an FT-IR spectrometer for a quasi in-situ chemical analysis using transmission mode: the gases are thus analyzed after cooling down to ambient temperature.

The operating conditions are typical of the CVI process [8]: (i) pressure between 1 and 10 kPa, (ii) temperature between 800 and 1300 K, (iii) total flow rate 100 to 500 sccm (std. cm 3 per second), (iv) initial dilution ratio [BCl 3 ]/[H 2 ] = 1/2 to 1/10. At each parameter set, data were collected at different times to check that a steady state is reached. Two independent mechanisms for relaxation towards steady state were identified: first, gaseous chemical equilibrium is reached a few minutes after the start of experiment; then, reactor thermal and mechanical steady-state achievement may take up to 20 minutes, and depends on the chosen temperature program. Not taking this last point into account may be a large source of error.

Numerical methods:

The modeling strategy starts with quantum chemical computations. First, equilibrium geometries, harmonic vibration frequencies and potential energy surfaces of a large set of chemical species in the system B/Cl/H are computed under the Density Functional Theory (DFT) framework at the B3Lyp/6-31G(d,p) level of theory using the Gaussian 2003 package [9]. This set includes H 2 , HCl, Cl 2 and BH i Cl j (0 ≤ i+j ≤ 4) species.

Vibration frequencies, known for being overestimated by DFT calculations, are scaled by the usual factor of 0.96. Partition functions, entropies and heat capacities are obtained from the well-known statistical equations [START_REF] Jensen | Introduction to Computational Chemistry[END_REF].

Central to the prediction of reaction rate constants is the calculation of the activation energy of the reaction, namely the potential energy difference between the transition state and the reactants (e.g. a variation of a few kcal/mol in the activation energy can change a reaction rate by several orders of magnitude). Calculation of energies thus requires numerical methods of high accuracy. In this work, we use the G3B3 method [START_REF] Curtiss | [END_REF], in which the energy is obtained from single point calculations (at the equilibrium B3Lyp/6-31G(d) geometry) at four different levels of theory (QCISD(T,E4T)/6-31G(d), MP4/6-31+G(d), MP4/6-31G(2df,p), MP2=Full/GTLarge), taking into account the electronic correlation to a large extent and a spin-orbit correction. Calculation with this method of ΔH f 0 (298K) of 155 molecules has given a root mean square deviation of only 1.37 kcal/mol [12].

Thermochemical computations:

The obtained heats of formations for all the investigated compounds are presented table 1.

Such a calculation requires the knowledge of atomic heats of formation. We use the values from JANAF tables [START_REF] Chase | JANAF Thermodynamic Tables, 3 rd ed[END_REF] for chlorine (28.98 kcal/mol) and hydrogen (52.07 kcal/mol) and the value recommended by Gurvich et al. [START_REF] Gurvich | Thermodynamic Properties of Individual Substances[END_REF] (134.93 kcal/mol) for boron. Experimental values (from [START_REF] Chase | JANAF Thermodynamic Tables, 3 rd ed[END_REF] and [START_REF] Gurvich | Thermodynamic Properties of Individual Substances[END_REF]) are given for comparison; the average standard deviation of the com-puted values with respect to these data shows that the agreement is excellent, especially with the Gurvich data. Discrepancies with JANAF tables concern primarily BH 2 , BCl 2 and BCl, most probably due to inaccuracies in these tables; indeed, similar computational results obtained with other high level quantum chemistry methods [12,[START_REF] Schlegel | [END_REF]16] ) in our model. We observed that all these species are in negligible quantities whatever the temperature.

Kinetic mechanism:

The determination of a kinetic mechanism first requires the calculation of the energy barriers associated to the chemical reactions which can take place in the gas phase. Looking at the typical CVD or CVI processing conditions, for T < 1500 K, we have concluded from the thermodynamic study that only a few light molecular species occur. As a consequence the reactions listed in table 2 are the most important candidates for the low temperature mechanism, furthermore considering that the energy barriers for their conversion are moderate (about 40 kcal/mol [7]).

We have thus reexamined Lengyel's mechanism [7] at the G3B3 level of theory (in their original work the DFT B3LYP/6-31G(2d,p) level was used). Transition states, energy barriers and partition functions have been characterized and the corresponding rate constants ap-proximated, according to Transition State Theory (TST). These reaction rates are fitted to an Arrhenius law; activation energies and pre-exponential factors are shown in table 2.

CVD reactor model :

The experimental reactor was first modeled to characterize the temperature profile along its axis. Total mass, momentum, and energy equations (featuring radiative heat transfer) were solved using a finite-volume software.

Figure 3 shows a typical temperature field obtained with 200 sccm of pure H 2 at 1 kPa and a susceptor temperature of 1300K. The radial temperature gradient is low and a onedimensional simulation represents a good approximation. Doubling the mass flow rate did not have a large influence: only a 2 mm lag on the axial temperature profile has been obtained. Also, composition did not have a large influence at high dilution because in every case H 2 is mainly responsible for heat conduction. For susceptor temperatures different from 1300 K it has been sufficient to consider a scaled profile based on the normalized 1300 K result. Once the temperature distribution is known, the species mass balance equations are solved using 1D-model and in steady state.

Results:

Figure 4 presents the steady state concentration of BCl 3 , BHCl 2 , HCl, BH 2 Cl and BH 3 at the reactor outlet obtained by the 1D model implementing our new kinetic mechanism. We see on Figure 4 that BCl 3 starts to decompose in a pronounced manner when the temperature exceeds 750 K. The main products of this decomposition are BHCl 2 and HCl. Around 800 K, BH 2 Cl starts to be observed in a non negligible amount. However, the ratio between BH 2 Cl and BHCl 2 concentrations remains around 1/100 for temperatures between 800 and 1400 K. Let us now compare the results obtained with the new mechanism to our experimental data, to results obtained by the mechanisms of Harris et al. [6] and Lengyel et al. [7] and to the thermodynamic equilibrium situation. We were not able to detect experimentally the presence of BH 2 Cl, which is coherent with Figure 4, showing that its concentration never exceeds 0.04 mol/cm we do not observe any additional IR signal we attribute this loss of boron to the beginning of the boron deposition process, in agreement with the observations of Sezgi et.al. [2]. Comparing concentrations obtained with our new model to experimental measurements for temperatures below the deposition limit, we see that the model agrees qualitatively well with experiments but that a shift of roughly 80 K is observed (the concentration given by the model at T equals the experimental concentration at (T-80) K). Of course our model does not predict the deposition process at higher temperatures as it is a purely homogeneous model. When the model of Lengyel et al [7] is used, similar trends can be observed but in that case the shift with experiments is around 180 K. This shows that increasing the quality of the ab-initio description permits to reduce the gap between experiments and calculations.

When the model by Harris et al. [6] is used, we see from figure 5 (respectively figure 6) that the decomposition of BCl 3 (respectively the apparition of BHCl 2 ) starts around 1200 K, that is, more than 400 K hotter than the experimental observations. Such a discrepancy can only be attributed to a mechanism which does not take into account the chemistry responsible for the decomposition of BCl 3 .

Concluding remarks:

CVD from BCl 3 and H 2 as precursors involves only three homogeneous molecular reactions at moderate temperature (T<1500K). Because of the relatively high activation energies of the radical initiation, the radical-containing route is negligible. Transition state theory was invoked at the G3B3 high level of theory to obtain reaction rates, and the resulting kinetic mechanism incorporated in a 1D model of our CVD reactor.

Theoretical steady state concentrations at the outlet of reactor agree well with experimental FTIR data at low temperature (T < 1000 K); the calculated values overestimate by only 80 K the evolution of the gas phase decomposition. This shift is about 180 K when using the mechanism of Lengyel et al. [7], which considers the same reaction but with data obtained at the lower DFT B3Lyp level of theory. Results obtained from the mechanism of Harris et al.

[6] dramatically differ from experimental data because it does not take into account the molecular reactions.

At higher temperatures both BCl 3 and BHCl 2 steady state concentrations are overestimated by our model. This is linked to an experimental loss of homogeneous boron species, which agrees well with the start of deposition process at about 1000 K.

This work is the starting point for a broader study concerning the deposition of boroncontaining carbides. It has been demonstrated here that the methodology is sound, as shown by the experimental validation; there remains to extend it to a larger gas-phase mechanism featuring carbon-containing species -for which validated kinetic databases already exist [17][18][19] -and to complete it with heterogeneous kinetic data, obtained either by fitting experimental deposition results to gas-phase concentrations, as in [20], or by quantum dynamic studies.

BH 3

 3 can only be observed at extremely low concentrations. It is hard to say whether the observed concentrations of BH 3 are real or due to numerical inaccuracy. Also displayed on Figure 5 are the concentrations of BCl 3 and BHCl 2 when only the first reaction of the kinetic mechanism (BCl 3 + H 2 -> BHCl 2 +HCl) is taken into account. We see that these concentrations are almost unaffected by the suppression of the two other reactions and that this single reaction is almost sufficient to properly describe BCl 3 decomposition in the investigated temperature range.

Figure 1 :

 1 Figure 1: Equilibrium molar fractions of molecular species in the gas phase as a function of temperature.

Figure 2 :

 2 Figure 2 : Equilibrium molar fractions of radicals in the gas phase as a function of temperature.

Figure 3 :

 3 Figure 3: Temperature field obtained with 200 sccm pure H 2 at 1 kPa and a susceptor temperature of 1300 K.

Figure 4 :

 4 Figure 4: Steady state concentrations at the reactor outlet as a function of temperature obtained from a 1D model of a CVD furnace using our new kinetic mechanism, either with a 3reaction model (lines), or with a 1-reaction model (symbols).

Figure 5 :

 5 Figure 5: Steady state BCl 3 concentration at the reactor outlet as a function of the susceptor temperature.

Figure 6 :

 6 Figure 6: Steady state BHCl 2 concentration at the reactor outlet as a function of the susceptor temperature.

  3 . Results for BCl 3 and BHCl 2 are presented in figures 5 and 6 respectively. Experimental concentrations of BCl 3 and BHCl 2 are related to the IR signals with the hypothesis of complete transformation of BCl 3 in BHCl 2 . It is seen first that the experimental decomposition of BCl 3 around 750 K coincides with the apparition of BHCl 2 situated between 600 and 800 K

(the lowest temperature at which we observed BHCl 2 is 773 K). We also see on those figures that around 1000 K both BCl 3 and BHCl 2 concentrations are decreasing with temperature. As

Table 1 :

 1 Room temperature heats of formation Δ f H 298K (kcal/mol) of BCl i H j (1 ≤ i+j ≤ 3) calculated at the G3B3 level of theory together with experimental values taken from JANAF tables a[START_REF] Chase | JANAF Thermodynamic Tables, 3 rd ed[END_REF] and from Gurvich et al. b[START_REF] Gurvich | Thermodynamic Properties of Individual Substances[END_REF] (available errors are given between parentheses). Absolute difference between theoretical and experimental value are indicated in Δ columns (kcal/mol).

	Δ f H 298K G3B3 Experiments a	Experiments b
		Values	|Δ|	Values	|Δ|
	BCl 3	-97.58 -96.3(0.5) 1.28 -96. 68(0.31) 0.9
	BHCl 2	-60.36 -59.3(1) 1.06 -60.2	0.16
	BH 2 Cl	-19.6		-19.3(4.8)	0.3
	BH 3	23.17 25.5(2.4) 2.33 21(2.4)	2.17
	BCl 2	-7.55 -19(3)	11.45
	BClH	32.59 27.5	5.09 31(5)	1.59
	BH 2	76.19 48(15)	28.19 76.1(2.6)	0.09
	BCl	41.58 33.8(4)	7.78 41(6)	0.58
	BH	103.48 105.8(2) 2.32 106.6(1.7)	3.12
	|Δ| Average		7.44	1.11

Table 2 :

 2 Activation energy E a (kcal/mol) and pre-exponential factor A (10 13 cm 3 /mol/s) for reactions(1)(2)(3). Values from Lengyel et al.[7] are indicated in parentheses for comparison.

	E a	A

Acknowledgements:

The Gaussian03 computations have been carried out at the M3PEC massively parallel computing facility of University Bordeaux 1. G. Reinisch acknowledges Snecma Propulsion Solide (Safran group) for PhD funding and wishes to thank R. Méreau (LPCM, University Bordeaux 1), C. Descamps (Snecma Propulsion Solide) for helpful discussions and S. Léger for help in CFD computations.


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

Bernard Vau
email: bernard.vau@ixblue.com

Tudor-Bogdan Airimitoaie
email: tudor-bogdan.airimitoaie@u-bordeaux.fr

Gabriel Buche
email: gabriel.buche@gipsa-lab.grenoble-inp.fr

Improving performance of adaptive feedforward noise attenuators using a dynamic adaptation gain

Keywords: active noise control, adaptive feedforward compensation, Youla-Kučera parametrization, adaptation algorithms, learning algorithms

The paper explores in detail the use of dynamic adaptation gain/learning rate (DAG) for improving the performance of adaptive feedforward attenuation schemes. The DAG is an ARMA (poles-zeros) filter embedded in gradient type adaptation/learning algorithms and generalizes the various improved gradient algorithms available in the literature. After introducing the DAG algorithm in the context of adaptive feedforward attenuation schemes and providing relationships with other algorithms, its design is developed. Strictly Positive Real (SPR) conditions play an important role in the design of the DAG. Then the stability issues for adaptive/learning systems using a DAG are discussed. The potential of the DAG is then illustrated by experimental results obtained on a relevant adaptive active noise control system.

Introduction

Many parameter adaptation algorithms (PAA) used in active noise and vibration control (ANVC) are based on the gradient technique [START_REF] Elliott | Signal processing for active control[END_REF]. There are a number of algorithms resulting either from the choice of the criterion to be minimized or from the approximation of the gradient. Other adaptation algorithms used in ANVC result from using a stability point of view for their development [START_REF] Landau | Adaptive and Robust Active Vibration Control-Methodology and Tests[END_REF]. Some of these stability-based algorithms can be interpreted as gradient type algorithms.

One of the most common features in applications is the use of a constant adaptation gain/learning rate in order to keep, on one hand, the alertness of the algorithm in the presence of changes in the characteristics of the disturbances (noise or vibration) to be attenuated, and on the other hand, to have a low computational load. Stability issues not addressed in many cases often lead to the use of very low values for the adaptation gains, resulting in very slow adaptation transients and unsatisfactory performance. The choice of the adaptation gain value is a major problem for applications.

Stability issues have been well understood in adaptive feedforward ANVC schemes. Appropriate algorithms which guarantee the stability of the system even for high adaptation gain values have been developed [START_REF] Landau | Adaptive feedforward compensation algorithms for active vibration control with mechanical coupling[END_REF][START_REF] Landau | IIR Youla-Kučera parameterized adaptive feedforward compensators for active vibration control with mechanical coupling[END_REF]. These algorithms on one hand take into account the inherent internal positive feedback occurring in adaptive feedforward attenuation schemes and on the other hand use the "a posteriori" residual error as adaptation error (leading to the normalization of the measured residual error). The implementation of these algorithms requires identification of the secondary and reverse paths for design of the appropriate filters acting on the measurements.

Once the stability issue is resolved, a significant problem remains to be addressed: How to choose the value of the adaptation gain assuring the best compromise between the speed of adaptation and steady state performance? More specifically, given a value of the adaptation gain assuring good steady state performance, how can one improve the adaptation speed? The answer to this question is one of the main objectives of this paper.

Solutions derived from the recursive least squares algorithm with forgetting factor have been proposed ( [START_REF] Landau | Adaptive and Robust Active Vibration Control-Methodology and Tests[END_REF][START_REF] Huang | A variable step-size fxlms algorithm for narrowband active noise control[END_REF]). However this approach leads to a significant increase in computational load ( [START_REF] Huang | A variable step-size fxlms algorithm for narrowband active noise control[END_REF]). One of these solutions is the so-called "constant trace algorithm" [START_REF] Landau | Adaptive and Robust Active Vibration Control-Methodology and Tests[END_REF][START_REF] Landau | Adaptive control, 2nd Edition[END_REF] implementing a time varying matrix adaptation gain 1 . Experimental results in ANVC show a certain improvement in adaptation transients but the computational load becomes much more important than in the case of constant adaptation gains. The signal processing community has also proposed algorithms using a time varying matrix adaptation gain as well as different simplifications in order to try to reduce the computational load (variable step size algorithms) [START_REF] Huang | A variable step-size fxlms algorithm for narrowband active noise control[END_REF]. The question is: does another solution for improving the adaptation transients exist that does not lead to a significant increase in computational load?

Recently, the concept of "dynamic adaptation gain" (DAG) has been introduced in [START_REF] Landau | Does a general structure exist for adaptation/learning algorithms?[END_REF] and simulation results (based on an academic example) indicate a high potentiality of this approach. The basic idea behind DAG can be summarized as follows: the gradient (or its approximation) is first filtered before it corrects the current estimated parameters through the (stationary) adaptation gain. It results from [START_REF] Landau | Does a general structure exist for adaptation/learning algorithms?[END_REF], that a filter with 3 coefficients allows significant improvement of the adaptation transients (the ARIMA2 algorithm in [START_REF] Landau | Does a general structure exist for adaptation/learning algorithms?[END_REF]). This dynamic adaptation gain can be interpreted also as a PID (with a filtered derivative) filter acting on the gradient2 .

The objectives of this paper are:

• To investigate the potentiality of dynamic adaptation gain in adaptive feedforward noise attenuation schemes using constant adaptation gain/learning rate.

• To provide experimental results obtained on a relevant test-bench dedicated to active noise control (ANC).

• To provide a methodology for design of dynamic adaptation gains.

As it will be shown in the paper, "passivity" concepts will play a key role in the design of dynamic adaptation gains.

The paper is organized as follows: the experimental set-up is described in Section 2. Then a basic adaptive feedforward configuration using a gradient type PAA is presented in Section 3. The dynamic adaptation gain/learning rate is introduced in the context of adaptive feedforward compensators in Section 4. The design of the dynamic adaptation gain is discussed in Sections 5 and 6. Experimental results obtained on the ANC test-bench introduced in Section 2 are presented in Section 7.

Experimental setup

The view of the test-bench used for experiments and its detailed diagram are shown in Fig. 1. The actual dimensions of the test-bench are given in Fig. 2.

The speaker used as the source of disturbances is labelled 1, while the control speaker is marked 2. At the pipe's open end, the microphone that measures the system's output transfer function between microphones (4→3) is called Primary Path (D), while the transfer function between the control speaker and the microphone (2→3) is denoted Secondary Path (G). The internal coupling found between (2→4) is denoted Reverse Path (M). These marked paths have a double differentiator behavior, since as input we have the voice coil displacement and as output the air acoustical pressure.

In this configuration, speakers are isolated inside wood boxes filled with special foam in order to create anechoic chambers and reduce the produced radiation noise. These boxes have dimensions 0.15 m × 0.15 m × 0.12 m, giving a chamber volume of 2.7 L.

Speakers and microphones are connected to a target computer with Simulink Real-time ® . A second computer is used for development and operation with Matlab. The sampling frequency is f s = 2500 Hz.

The various paths are described by models of the form3 :

X(q -1 ) = q -dx B X (q -1 ) A X (q -1 ) = q -dx b X 1 q -1 + ... + b X n B X q -n B X 1 + a X 1 q -1 + ... + a X n A X q -n A X (1)
with

B X = q -1 B * X for any X ∈ {G, M, D}. Ĝ = q -d G BG A G , M = q -d M BM A M , and D = q -d T BT A T
denote the identified (estimated) models of G, M , and D. The system's order is defined as (the indexes G, M , and T have been omitted):

n = max(n A , n B + d).
The frequency characteristics of the identified models for the primary4 , secondary and reverse paths are shown in Fig. 3. These characteristics present multiple resonances (low damped complex poles) 5 and anti-resonances (low damped complex zeros). One can see that the secondary path has a high gain between 70 to 250 Hz, which means that disturbances can be efficiently attenuated in this region. It is also clear that the reverse path has a significant gain on a large frequency range (up to 400 Hz) so its effect cannot be neglected. These models have been identified experimentally using the techniques described in [START_REF] Melendez | Data driven design of tonal noise feedback cancellers[END_REF]. The orders and the pure delays of the various identified models are given in Table 1.

It is important to note that most of the implementations of the adaptive feedforward compensation systems are close to a co-location of the residual noise measurement and of the secondary source used for compensation. See for example [START_REF] Zeng | Recursive filter estimation for feedforward noise cancellation with acoustic coupling[END_REF], [START_REF] Fraanje | Robust and fast schemes in broadband active noise and vibration control[END_REF], [START_REF] Hu | Feedforward active noise controller design in ducts without independent noise source measurements[END_REF]. A ratio of 1/6 to 2/6 between the length of the secondary path and the length of the primary path is Model [START_REF] Gaudio | Connections between adaptive control and optimization in machine learning[END_REF] 19 0 Secondary (G) [START_REF] Anderson | Stability of adaptive systems[END_REF] 24 5 Reverse (M) 22 25 5

n B n A d Primary (D)
Table 1: Orders of the identified system paths.

used in these references. Nevertheless, there are new potential applications areas where technological constraints will not allow a configuration to be close to a co-location. For this reason the ratio between the length of secondary path and the length of the primary path has been chosen slightly larger than 4/6. The objective of the active control system is to attenuate incoming unknown broad-band noise disturbances.

The basic adaptive feedforward configuration

For the purpose of this paper (comparison of various adaptation/learning algorithms), the FIR Youla-Kučera configuration for the feedforward compensator has been selected. This configuration has been described previously in the active vibration context [START_REF] Landau | A Youla-Kučera parametrized adaptive feedforward compensator for active vibration control with mechanical coupling[END_REF], [START_REF] Landau | IIR Youla-Kučera parameterized adaptive feedforward compensators for active vibration control with mechanical coupling[END_REF] and in the active noise context [START_REF] Landau | Robust and adaptive feedback noise attenuation in ducts[END_REF][START_REF] Airimitoaie | Algorithms for Adaptive Feedforward Noise Attenuation-A Unified Approach and Experimental Evaluation[END_REF]. It offers the great advantage of defining a priori the closed loop poles of the internal positive feedback loop independently of the values of the adaptive FIR filter assuring the attenuation of the incoming disturbance. For a detailed description of the algorithm and of the stability conditions, see for example [START_REF] Landau | Robust and adaptive feedback noise attenuation in ducts[END_REF]. In these papers, a stability approach has been considered for the synthesis of the PAAs. In this paper, the PAAs are synthesized using the "gradient rule"6 and then the stability of the system is analyzed. This approach is necessary for the subsequent introduction of the new algorithms which can be viewed as an improvement of the "gradient" type algorithms.

The corresponding block diagram for the adaptive feedforward noise compensation using the FIR Youla-Kučera (FIR-YK) parametrization of the feedforward compensator is shown in Figure 4. ŷ(t) denotes the effective output of the measurement device (microphone in ANC) providing an image of the incoming disturbance and which will serve as input to the adaptive feedforward filter. The output of this filter, denoted by û(t), is applied to the actuator through an amplifier. The transfer function G (the secondary path) characterizes the dynamics from the output of the feedforward filter to the residual noise measurement (amplifier + actuator + dynamics of the system). The transfer function D characterizes the primary path. The transfer function M characterizes the positive feedback path.

The adjustable filter Q has the structure: and the parameters q i will be adapted in order to minimize the residual error. One defines

Q(q -1 ) = q0 + q1 q -1 + ... + qn Q q -n Q (2)
w T = [q 0 , q 1 , q 2 , . . . , q n Q ] (3) ŵT = [q 0 , q1 , q2 , . . . , qn Q ] (4) u T (t) = [α(t + 1), α(t), . . . , α(t -n Q + 1)] (5)
where:

α(t + 1) = B M û(t + 1) -A M ŷ(t + 1) = B * M ŷ(t) -A M û(t + 1) (6)
The a priori output of the feedforward compensator is given by: û• (t + 1) = -S * 0 û(t) + R 0 ŷ(t + 1) + ŵT (t)u(t) [START_REF] Landau | Does a general structure exist for adaptation/learning algorithms?[END_REF] where û(t), û(t -1), ... are the "a posteriori" outputs of the feedforward filter generated by

û(t + 1) = -S * 0 û(t) + R 0 ŷ(t + 1) + ŵT (t + 1)u(t) (8)
The measured input to the feedforward filter satisfies the following equation (when feedforward compensation is active)

ŷ(t + 1) = d(t + 1) + B * M (q -1 ) A M (q -1 ) û(t). ( 9
)
The unmeasurable value of the output of the primary path is denoted x(t). The unmeasurable a priori output of the secondary path will be denoted ẑ• (t + 1).

ẑ• (t + 1) = ẑ(t + 1| ŵ(t)) = B * G (q -1 ) A G (q -1 ) ŷ(t) (10)
The a posteriori unmeasurable value of the output of the secondary path is denoted:

ẑ(t + 1) = ẑ(t + 1| ŵ(t + 1)) (11)
The a priori adaptation error is defined as:

ν • (t + 1) = ν(t + 1| ŵ(t)) = -e • (t + 1) = -x(t + 1) -ẑ• (t + 1) (12)
where e • (t+1) is the measured residual noise. The a posteriori adaptation error (computed) will be given by:

ν(t + 1) = ν(t + 1| ŵ(t + 1)) = -x(t + 1) -ẑ(t + 1). ( 13
)
When using an estimated filter N with constant parameters: ŷ

• (t) = ŷ(t), ẑ• (t) = ẑ(t) and ν • (t) = ν(t).
Using the results of [START_REF] Landau | A Youla-Kučera parametrized adaptive feedforward compensator for active vibration control with mechanical coupling[END_REF], the a posteriori residual adaptation error can be expressed as:

ν(t + 1) = ν(t + 1| ŵ(t + 1)) = H(q -1 )[w -ŵ(t + 1)] T u(t) (14)
where7

H(q -1 ) = A M G P 0 ; P 0 = A M S 0 -B M R 0 (15)
Assume that we would like to minimize at each step a quadratic criterion in terms of the a posteriori residual error min ŵ(t+1)

J(t + 1) = [ν(t + 1)] 2 (16)
A solution can be provided by the gradient rule. The corresponding parameter adaptation/learning algorithm (PALA) will have the form:

ŵ(t + 1) = ŵ(t) -F ▽ w J(t + 1) = ŵ(t) -F ∂J(t + 1) ∂ ŵ(t + 1) ( 17
)
where F is the matrix adaptation gain/learning rate and ▽ w J(t + 1) is the partial gradient of the criterion given in [START_REF] Landau | Adaptive control : the model reference approach[END_REF] with respect to ŵ(t+1) which has the expression ∂J(t+1) ∂ ŵ(t+1) . There are two possible choices for the matrix adaptation gain/learning rate: (i) F = γI; γ > 0; (ii) F > 0 (positive definite matrix). For the remainder of the paper we will use the option F = γI. The term (constant) adaptation gain or learning rate is used for γ. Neglecting the non-commutativity of the transfer operator (valid for low values of the adaptation gain), ( 14) can be re-written as:

ν(t + 1) = ν(t + 1| ŵ(t + 1)) = [w -ŵ(t + 1)] T H(q -1 )u(t) (18)
Using this equation, the expression of the partial gradient becomes:

▽ w J(t + 1) = -H(q -1 )u(t)ν(t + 1) (19)
Unfortunately, H(q -1 ) is not known exactly and should be replaced by its estimation Ĥ(q -1 ) based on the identified models of the secondary and reverse paths. Therefore the gradient of the criterion will be computed using:

▽w J(t + 1) = L(q -1 )u(t)ν(t + 1); L(q -1 ) = Ĥ(q -1 ) (

One defines also the regressor vector (a filtered observation vector) as:

r(t) = L(q -1 )u(t) = [α f (t + 1), α f (t), . . . , α f (t -nQ + 1)] (21)
where:

α f (t + 1) = L(q -1 )α(t + 1) (22)
and the basic gradient type PALA will take the form:

ŵ(t + 1) = ŵ(t) + γr(t)ν(t + 1), (23)
Using R 0 = 0 and S 0 = 1, the poles of the internal positive closed loop will be defined by A M (see Eq. ( 15)) and they will remain unchanged (see [START_REF] Landau | Robust and adaptive feedback noise attenuation in ducts[END_REF] for details). In this context, the filter used in [START_REF] Anderson | Stability of adaptive systems[END_REF] will be L = Ĝ. Defining the parameter error:

w(t) = ŵ(t) -w (24)
one obtains from ( 14) and [START_REF] Bourlès | Linear Systems[END_REF]:

ν(t + 1) = - H(q -1 ) Ĥ(q -1 ) w(t + 1)] T r(t) (25)
w(t + 1) = w(t) + γr(t)ν(t + 1) (26)

w(t + 1) T r(t) = w(t) T r + γr T (t)r(t)ν(t + 1) (27)
Equations ( 25) and ( 27) define an equivalent feedback system shown in Fig. 5 and the associated linear transfer operator appearing in the equivalent feedforward path is:

H ′ (q -1 ) = H(q -1 ) L(q -1 ) (28)
Using the results of [6, Theorem 3.2] one concludes that asymptotic stability will be obtained for any values of the adaptation gain γ and any initial conditions if H ′ (z -1 ) is a strictly positive real (SPR) transfer function. This is a very mild condition as far as a good experimental identification of the models is carried out.

To implement the algorithm, we need an expression of ν(t + 1) as a function of ν • (t + 1). From [START_REF] Pouyanfar | A survey on deep learning: Algorithms, techniques, and applications[END_REF], one gets: where H(q -1 ) = 1 + q -1 H * (q -1 ); Ĥ(q -1 ) = (1 + q -1 Ĥ * (q -1 ) (30)

ν(t + 1) = [w -ŵ(t + 1)] T r(t) -Ĥ * (q -1 )ν(t) + H * (q -1 ) [w -ŵ(t)] T r(t -1) (29)
The a priori adaptation error obeys the equation

ν • (t + 1) = [w -ŵ(t)] T r(t) -H * (q -1 ) [w -ŵ(t)] T r(t -1) -Ĥ * (q -1 )ν(t) (31)
The a posteriori adaptation/learning error can then be written as:

ν(t + 1) = ν • (t + 1) -[ ŵ(t + 1) -ŵ(t)] T r(t) = ν • (t + 1) -γr(t) T r(t)ν(t + 1) (32)
which leads to:

ν(t + 1) = ν • (t + 1) 1 + γr T (t)r(t) (33)
For the purpose of this paper, it is convenient to express the adaptation algorithm as:

r(t) = u f (t) = L(q -1 )u(t) (34)
ν(t + 1) = ν • (t + 1) 1 + γr T (t)r(t) (35) ŵ(t + 1) = ŵ(t) + γr(t)ν(t + 1) (36) ( 37
)
where ν • (t + 1) is given in [START_REF] Landau | A Youla-Kučera parametrized adaptive feedforward compensator for active vibration control with mechanical coupling[END_REF]. The asymptotic stability condition for any values of the adaptation gain γ and any initial conditions ŵ(0), ν • (0) is that H ′ (z -1 ) given in (28) should be strictly positive real (SPR).

An adaptation algorithm using a dynamic adaptation gain/learning rate

The new adaptation algorithm which will be discussed in this section will have the form8

ŵ(t + 1) = ŵ(t) + γ C(q -1 ) D ′ (q -1 ) [r(t)ν(t + 1)] (38)
where C(q -1 ) D ′ (q -1 ) is termed the "dynamic adaptation gain/learning rate" (DAG). It is in fact a MIMO diagonal transfer operator having identical terms. All the diagonal terms are of the form:

H ii DAG (q -1 ) = C(q -1 ) D ′ (q -1 ) = 1 + c 1 q -1 + c 2 q -2 + .. + c n C q -n C 1 -d ′ 1 q -1 -d ′ 2 q -2 -..d ′ n D ′ q -n ′ D ( 39
)
The effective implementation of the algorithm given in (38) leads to:

ŵ(t + 1) = d 1 ŵ(t) + d 2 ŵ(t -1) + . . . + d n D ŵ(t -n D ) + γ[r(t)ν(t + 1) + c 1 r(t -1)ν(t) + c 2 r(t -2)ν(t -1) + . . . + c n C r(t -n C )ν(t -n C + 1)] ( 40
)
where:

d i = (d ′ i -d ′ i-1 ) ; i = 1, ...n D ; d ′ 0 = -1, d ′ n D = 0 (41)
To implement the algorithm one needs a computable expression for ν(t + 1). Using a similar development as for the basic gradient adaptation algorithm one finds:

ν(t + 1) = ν • (t + 1) 1 + γr T (t)r(t) (42)
Relations with other algorithms Many algorithms have been proposed for improving the "gradient rule" in order to accelerate the convergence of the adaptation/learning algorithm for a given adaption gain/learning rate. The algorithm of (38) is termed ARIMA (Autoregressive with Integrator Moving Average ) algorithm [START_REF] Landau | Does a general structure exist for adaptation/learning algorithms?[END_REF]. The various algorithms described in the literature are of MAI form or ARI form. The MAI form includes "Integral + Proportional" algorithm [START_REF] Landau | Adaptive control, 2nd Edition[END_REF][START_REF] Airimitoaie | Improving adaptive feedforward vibration compensation by using integral+proportional adaptation[END_REF] (c 1 ̸ = 0, c i = 0, ∀ i > 1, d ′ i = 0, ∀ i ≥ 0), "Integral + Proportional + Derivative" algorithm [START_REF] Landau | Adaptive control : the model reference approach[END_REF][START_REF] Landau | Does a general structure exist for adaptation/learning algorithms?[END_REF] [START_REF] Schmidt | Stochastic average gradient[END_REF][START_REF] Pouyanfar | A survey on deep learning: Algorithms, techniques, and applications[END_REF]. The ARI form includes "Conjugate gradient" and "Nesterov" algorithms [START_REF] Livieris | A survey on algorithms for training artificial neural networks[END_REF][START_REF] Gaudio | Connections between adaptive control and optimization in machine learning[END_REF] (c i = 0, f or i = 1, 2, .., d ′ 1 ̸ = 0, d ′ i = 0, f or i > 1) as well as the "Momentum back propagation" algorithm [START_REF] Jacobs | Increased rates of convergence through learning rate adaptation[END_REF] which corresponds to the conjugate gradient plus a normalization of γ by (1 -d ′ 1 )9 . For more details see [START_REF] Landau | Does a general structure exist for adaptation/learning algorithms?[END_REF]. Application of the I+P algorithm in active vibration control is discussed in [START_REF] Airimitoaie | Improving adaptive feedforward vibration compensation by using integral+proportional adaptation[END_REF].

(c 1 , c 2 , c i = 0 f or i > 2, d ′ i = 0, i ≥ 0), "Averaged gradient" (c i , f or i = 1, 2, ..., d ′ i = 0, ∀ i ≥ 0)
A particular form of the ARIMA algorithms termed "ARIMA2" (c

1 , c 2 , c i = 0, ∀ i > 2, d ′ 1 ̸ = 0, d ′ i = 0, ∀ i > 1)
will be studied subsequently and evaluated experimentally. 10

Design of the dynamic adaptation gain -performance issues

The dynamic adaptation gain/learning rate will introduce a phase distortion on the gradient depending on the frequency.

• Assume that the algorithms should operate for all frequencies in the range: 0 to 0.5f s .

• Assume that the gradient of the criterion to be minimized contains a single frequency.

In order to minimize the criterion, the phase distortion introduced by the dynamic adaptation gain/learning rate should be less than 90 • at all the frequencies. In other terms, the transfer function

H DAG = C(z -1 )
D ′ (z -1 ) should be SPR. But a transfer function which is SPR has its poles and zeros inside the unit circle (asymptotically stable poles and zeros). For transfer functions having their poles and zeros inside the unit circle one has the following property:

DAG Property 1: Assume that the polynomials C(z -1 ) and D ′ (z -1 ) have all their zeros inside the unit circle, then:

I = π 0 log C(e -iω ) D ′ (e -iω ) dω = 0 (43)
The proof of this property is given in Appendix A. This result allows to conclude that the average gain (in dB) over the frequency range 0 to 0.5f s is 0 dB (=1), i.e. on average this filter will not modify the adaptation gain/learning rate. It is just a frequency weighting of the adaptation gain/learning rate. To be more specific, Figure 6 shows the frequency characteristics of two DAGs which will be subsequently used in the experimental section 11 . It can be observed first that the phase is within the range ±90 • , i.e. they are SPR. Then one can observe that effectively the average gain over the frequency range 0 to 0.5f s (f s = 2500 Hz) is 0 dB. Now examining the magnitude, one observes that both are low pass filters amplifying low frequencies. This means that if the frequency content of the gradient is in the low frequency range, the effective adaptation gain/learning rate will be larger than γ, which should have a positive effect upon the adaptation/learning transient. In particular the DAG which has a larger gain in low frequencies (ARIMA2) should provide better performance than the (I+P) DAG (which is indeed the case-see Section 7). Furthermore, since the average gain is 0 dB, the filters shown in Figure 6 will introduce an attenuation in the high frequency range, reducing the effective adaptation gain in this region. This will reduce the influence of the measurement noise located usually in the high frequency region and steady state performance is improved (in general) with respect to the gradient (see Section 7).

Since we need to have a DAG which is SPR, we will provide subsequently the tools for designing a SPR DAG. We will consider the case of the ARIMA2 algorithm introduced in [START_REF] Landau | Does a general structure exist for adaptation/learning algorithms?[END_REF]. The DAG in this case will have the form:

H DAG = C(q -1 ) D ′ (q -1 ) = 1 + c 1 q -1 + c 2 q -2 1 -d ′ 1 q -1 (44)
A criterion for the selection of c 1 , c 2 and d ′ 1 in order that the DAG be SPR is given below.

SPR criterion for DAG:

The conditions assuring that

H DAG (z) = 1+c 1 z -1 +c 2 z -2 1-d ′ 1 z -1
is strictly positive real (SPR) are: • for c 2 ≤ 0, c 1 must be such that

-1 -c 2 < c 1 < 1 + c 2 • for c 2 ≥ 0 -if the following condition is satisfied 2(d ′ 1 -c 2 ) < 2(c 2 -c 2 2 )(1 -d ′ 2 1 ) < 2(d ′ 1 + c 2 )
the maximum bound on c 1 is given by

c 1 < d ′ 1 -3d ′ 1 c 2 + 2 2(c 2 -c 2 2 )(1 -d ′ 2 1 )
otherwise the maximum bound on c 1 is given by

c 1 < 1 + c 2 -if the following condition is satisfied 2(d ′ 1 -c 2 ) < -2(c 2 -c 2 2 )(1 -d ′ 2 1 ) < 2(d ′ 1 + c 2 )
the minimum bound on c 1 is given by

c 1 > d ′ 1 -3d ′ 1 c 2 -2(c 2 -c 2 2 )(1 -d ′ 2 1 )
otherwise the minimum bound on c 1 is given by

c 1 > -1 -c 2
The proof of this result is given in Appendix B. From these conditions, closed contours in the plane c 2 -c 1 can be defined for the different values of d ′ 1 allowing to select c 1 and c 2 for a given value of d ′ 1 so that the DAG is SPR 12 . Such a diagram is presented in Figure 7.

Design of the dynamic adaptation gain-stability issues

Eq. ( 38) can be expressed also as

ŵ(t + 1) = H P AA (q -1 )γr(t)ν(t + 1) (45)
where H P AA is a MIMO diagonal transfer operator having identical terms. All the diagonal terms are described by:

H ii P AA (q -1 ) = 1 + c 1 q -1 + c 2 q -2 + .. + c n C q -n C (1 -q -1 )(1 -d ′ 1 q -1 -d ′ 2 q -2 -..d ′ n D ′ q -n ′ D ) = C(q -1 ) (1 -q -1 )D ′ (q -1 ) = C(q -1 ) D(q -1 ) (46)
The relation between the coefficients of polynomials D and D' is given in (41). One has the following result for the stability of the adaptive feedforward attenuation scheme using a dynamic adaptation gain/learning rate: Stability Conditions with DAG: For the system described by Eqs (2) through (14) using the PALA of (40), ( 41) and (42), one has lim t→∞ ν(t + 1) = 0 for any positive adaptation gain γ and initial conditions ŵ(0), ν(0), if:

• The transfer function associated with the transfer operator H ii P AA (q -1 ) given in (46) is a positive real (PR) transfer function with a pole at z = 1.

• The transfer function H ′ (z -1 ) given in (28) is a strictly positive real (SPR) transfer function The proof of this result is given in Appendix C.

For the particular case of the ARIMA2 algorithm, in order to operate with high adaptation/learning rate, the coefficients c 1 , c 2 and d ′ 1 should be chosen such that the DAG is SPR and the H ii P AA is PR, i.e.

H ii P AA = 1 + c 1 q -1 + c 2 q -2 1 -d 1 q -1 -d 2 q -2 = 1 + c 1 q -1 + c 2 q -2 (1 -q -1 )(1 -d ′ 1 q -1 ) (47)
should be PR. It is important for applications to give explicit bounds for the selection of the coefficients c 1 , c 2 , d ′ 1 in order to guarantee the positive realness of the embedded ARIMA filter. One has the following result:

PR conditions for H P AA : In order that the transfer operator H ii P AA given in (47) be characterized by a PR transfer function, the necessary and sufficient conditions are:

-1 < d ′ 1 < 1 (48) 0 ≤ δ ≤ 2 (49) -1 ≤ d ′ 1 - γ 1 -δ/2 ≤ 1 ( 50
)
where:

δ = 1 + c 1 + c 2 1 -d ′ 1 ; γ = d ′ 1 c 1 + d ′ 2 1 + c 2 d ′ 1 -1 (51)
The proof of this criterion is given in Appendix D.

From these conditions, closed contours in the plane c 2 -c 1 can be defined for the different values of d ′ 1 allowing to select c 1 and c 2 for a given value of d ′ 1 so that the H P AA is PR 13 . Such a diagram is presented in Figure 8. It is also interesting to see the intersections of the contours assuring the SPR of the H ii GAD with the contours assuring that H ii P AA is PR. Such an intersection is shown in Fig. 9. From this figure one can conclude that not all the SPR H GAD will lead to an H P AA PR. In such cases the performances are improved for low adaptation gain, but one cannot guarantee asymptotic stability for large values of the adaptation gain γ (see [START_REF] Landau | Does a general structure exist for adaptation/learning algorithms?[END_REF] for a simulated example). Fig. 9 shows also that there is a region where despite H P AA being PR, H DAG is not SPR. For such configurations, large adaptation gains can be used but the adaptation transient is slower than for the basic gradient algorithm. For small values of the adaptation gains/learning rates the passivity/stability condition can be relaxed using averaging [START_REF] Anderson | Stability of adaptive systems[END_REF]. Using the results of [START_REF] Landau | Adaptive feedforward compensation algorithms for active vibration control with mechanical coupling[END_REF], under the hypothesis of an input signal spanning all the frequencies up to half of the sampling frequency, passivity in the average will be assured if the frequency interval where H ii P AA is not positive real is smaller than the frequency interval where H ii P AA is positive real. In fact what is important is that the H ii P AA is PR in the frequency region of operation (mainly defined by the spectrum of the input signals to the systems).

Experimental results

Two types of signals used as unknown disturbances acting on the system have been considered:

• A broad-band disturbance 70-170 Hz, to which two tonal disturbances located at 100 Hz and 140 Hz have been added

• A broad-band disturbance 70-220 Hz

The steady state and transient attenuation 14 will be evaluated for the various values of the parameters c 1 , c 2 and d ′ 1 given in Tables 2 and3. The system will operate in open loop during the first 15 s. For all the experiments the constant adaptation gain γ has the value 0.2 and the adjustable filter Q has 60 parameters. 10 shows the time response of the system as well as the evolution of the global attenuation when using the Gradient algorithm (top) and the ARIMA2 (2) algorithm (bottom) (last row of Table 2). One can observe that the ARIMA2 (2) algorithm converges faster than the Gradient algorithm and provides also better attenuation in steady state. Figure 11 shows for the same experiments a comparison of the PSD of the residual noise in the absence of the control and the PSD of the residual noise at 600 s under the control effect. One can observe that the ARIMA2 algorithm provides better attenuation in the low frequencies as well as in the high frequency range when compared with the Gradient algorithm. This is coherent with the frequency characteristics of the ARIMA 2 dynamic adaptation gain, shown in Figure 12. It has a larger gain in low frequencies and introduces For |z| = 1, one has z = e jω , dz = je jω dω, and dz z = jdω. Thus one can write

I = 1 j T log 1 + nC k=1 c k z -k dz z - T log 1 - nD ′ k=1 d ′ k z -k dz z   (A.1)
where T is the unit circle.

On the other hand 1 + nC k=1 c k z -k = 1 + nC k=1 c k z k for |z| = 1. The poles of f (z) = log 1 + nC k=1 c k z k are the zeros of 1 + nC k=1 c k z k , and they all lie outside the unit 22 circle (by assumption, otherwise H P AA cannot be PR). Therefore the function f : z → f (z) is holomorphic in the open unit circle, and one can apply the Cauchy Integral's Formula (see [23, pg. 411]). This formula yields

1 j T f (z) dz z = Ind(T, 0)f (0)
where Ind(T, 0) is the index of the unit circle with respect to z = 0. One has Ind(T, 0) = 1, and f (0) = log(1) = 0. Therefore one gets

1 j T log 1 + nC k=1 c k z -k dz z = 1 j T f (z) dz z = 0
The same machinery can be applied mutatis mutandis for f (z) = log 1 Since the equivalent path is passive, Using Theorem C2 from [START_REF] Landau | Adaptive control, 2nd Edition[END_REF] one conclude that if the equivalent feedforward path is characterized by a SPR transfer function, one has that lim t→∞ ν(t + 1) = 0 for any initial conditions.

-nD ′ k=1 d ′ k z k

Appendix D. Proof of PR conditions for H P AA

A first necessary condition is given by condition (48) assuring that the poles of the transfer function H P AA are inside or on the unit circle. By performing a partial fraction expansion of H P AA , one has:

H P AA (q -1 ) = 1+ δq -1 1-q -1 + γq -1 1-d ′ 1 q -1
. Set β such that β ∈]0, 1[, one can write H P AA (q -1 ) = H 1 (q -1 ) + H 2 (q -1 ) with H 1 (q -1 ) = β + δq -1

1-q -1 = β 1-β-δ β q -1 1-q -1 and H 2 (q -1 ) = (1 -β)

+ γq -1 1-d ′ 1 q -1 = (1-β) 1-d ′ 1 -γ 1-β q -1 1-d ′ 1 q -1
. Since H 1 , H 2 are first order transfer function operators, and since β > 0, 1 -β > 0, a sufficient condition for H 1 and H 2 to be both PR is that their zeros be inside or on the unit circle. This is assured if the two following conditions are met simultaneously

-1 ≤ β -δ β ≤ 1 (a)

Figure 1 :

 1 Figure 1: Duct active noise control test-bench photo (top) and block diagram (bottom).

(

  residual noise e(t)) is denoted as 3. s(t) is the disturbance. Inside the pipe, close to the source of the disturbances, the second microphone, labelled 4, measures the disturbance's image, denoted as y(t). u(t) is the control signal. The transfer function between the disturbance speaker and the output microphone (1→3) is called Global Primary Path. The

Figure 2 :

 2 Figure 2: Duct active noise control test-bench dimension.

Figure 3 :

 3 Figure 3: Frequency characteristics of the Primary, Secondary and Reverse paths identified models.

Figure 4 :

 4 Figure 4: Feedforward AVC with FIR-YK adaptive feedforward compensator.

Figure 5 :

 5 Figure 5: Equivalent feedback system for the adaptive feedforward attenuation scheme.

Figure 6 :

 6 Figure 6: Frequency characteristics of two DAGs (used in the experiments).

Figure 7 :

 7 Figure 7: SPR contours in the plane c 2 -c 1 for various values of d ′ 1 (0, 0.5, 0.9) for the DAG given in Eq. (44).

Figure 8 :

 8 Figure 8: PR contours in the plane c 2 -c 1 for various values of d ′1 (0, 0.5, 0.9) for the H P AA of Eq. (47)

5 'Figure 9 :

 59 Figure 9: Intersection in the plane c 1 -c 2 of the contour H P AA = P R with the contour H GAD = SP R for d ′ 1 = 0.5

Figure 10 :

 10 Figure 10: Time evolution of the residual noise using the Gradient algorithm (top) and using the ARIMA2 algorithm (bottom).

Figure 11 :

 11 Figure 11: Comparison of the residual noise PSD in open loop and in closed loop (t=600 s) using the Gradient algorithm (top) and using the ARIMA2 (2) algorithm. Disturbance 70 -170 Hz + 2 tonal disturbances. (bottom).

Figure 12 :

 12 Figure 12: Frequency figcharacteristics of H DAG for the various algorithms given in table 3

Figure 13 :

 13 Figure 13: Time evolution of the global attenuation for the algorithms ofTable 2 (top) and of Table 3 (bottom). Disturbance 70 -170 Hz + 2 tonal disturbances.

Figure 14 :

 14 Figure 14: Time evolution of the global attenuation for the algorithms of Table 2 (top) and of Table 3 (bottom). Disturbance 70 -220 Hz.

Table 2 :

 2 ARIMA algorithms with SPR H DAG

	Algorithm	H P AA P R H DAG SP R	c 1	c 2	d ′ 1
	Gradient	Y	Y	0	0	0
	Conj.Gr/Nest.	N	Y	0	0	0.9
	I+P+D (γ P = -2γ D )	N	Y	1	0.8	0
	I+P	Y	Y	0.99	0	0
	ARIMA2 (1)	N	Y	-0.55 0.4 0.75
	ARIMA2 (2)	N	Y	0.99	0	0.9

Table 3 :

 3 ARIMA algorithms with SPR H DAG and H P AA PR.

	Algorithm	H P AA P R H DAG SP R c 1	c 2	d ′ 1
	Gradient	Y	Y	0	0	0
	Conj.Gr/Nest.	Y	Y	0	0	0.333
	I+P+D (γ P = -2γ D )	Y	Y	0.1 0.333	0
	I+P	Y	Y	0.99	0	0
	ARIMA 2	Y	Y	0.1 0.08 0.22
	7.1. Broad-band disturbance 70-170 Hz plus two tonal disturbances located at 100 Hz and
	140 Hz
	Figure

Table 2 (

 2 top) and of Table 3 (bottom). Disturbance 70 -220 Hz. and since the the function C(e -iω ) D ′ (e -iω ) is even one gets the claimed result. Appendix C. Stability of the system using the dynamic adaptation gain Consider Eq. (45). Using the [A, B, C, D] state space representation associated to H P AA (z) one gets: Eqs. (C.1), and (C.3) define an equivalent feedback system, the equivalent feedback path being defined by (C.1) and (C.3). In order to apply passivity theorem the equivalent feedback path should be passive. One can use Theorem 3.3.1 from [6]: Theorem 1. For a PAA having the form of Eqs (C.1) and (C.2), the equivalent feedback path described by Eqs. (C.1), and (C.3) is passive, i.e., η(0, t 1 ) =

	x(t + 1) = Ax(t) + Br(t)ϵ(t + 1)	(C.1)
	w(t + 1) = Cx(t) + Dr(t)ν(t + 1)	(C.2)
	and respectively:
	r T (t) w(t + 1) = r T (t)Cx(t) + r T (t)Dr(t)ν(t + 1)	(C.3)
	t 1

t=0 ν(t + 1)r T (t) w(t + 1) ≥ -γ 2 ; γ 2 < ∞ for all t ≥ 0 (C.

4) if the associated linear system [A, B, C, D] described by Eqs.(C.1) and (C.2) is passive, or equivalently, if H P AA (z) given in Eq. (45) is a PR transfer matrix.

The trace of the matrix adaptation gain is kept constant and it is equal to the desired adaptation gain per parameter multiplied by the number of parameters.

The correcting term is no more only proportional to the gradient but it depends also upon the integral of the gradient as well as upon the filtered variations of the gradient.

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.

The primary path model has been exclusively used for simulation purposes.

The lowest damping is around 0.01.

In the gradient rule the correction of the parameters is done in the direction of the gradient of the criterion to be minimized with the minus sign.

P 0 defines the poles of the internal positive closed loop formed by the feedforward compensator and the reverse path.

It can be interpreted as a dynamic gradient rule.

There are very few indications on how to choose the various weights in the above mentioned algorithms.

The algorithms mentioned above can be viewed as particular cases of the ARIMA2 algorithm.

ARIMA2 filter with c 1 = 0.99, c 2 = 0, d ′ 1 = 0.9 and I+P filter with c 1 = 0.99, c 2 = 0, d ′ 1 = 0.

A Matlab routine is available for drawing these contours.

A Matlab routine is available for drawing these contours.

The attenuation is defined as the ratio between the variance of the residual noise in the absence of the control and the variance of the residual noise in the presence of adaptive feedforward compensation. The variance is evaluated over a horizon of 3 s.

), otherwise this lower bound is given by c 1 > -c

-1. This ends the proof.

attenuation in the high frequencies (in other terms: the adaptation gain is higher in low frequencies and lower in the high frequencies).

The top of Fig. 13 shows the time evolution of the global attenuation for the algorithms considered in Table 2. As it can be seen, there is a clear improvement in the adaptation transient using ARIMA2 [START_REF] Landau | Adaptive and Robust Active Vibration Control-Methodology and Tests[END_REF] (last row of Table 2) with respect to the Gradient algorithm (first row of Table 2). The Conjugate Gradient/Nesterov algorithm also provides a significant improvement of the transient. The ARIMA2 (2) algorithm reduces the adaptation/learning transient by a factor of three and a half. One also sees an improvement in steady state attenuation with respect to gradient adaptation. The other algorithms (from Table 2) also provide an improvement with respect to the gradient algorithms. Note that except for the I+P algorithm, H P AA is not PR, therefore there is a risk of instability for large values of the adaptation gain γ.

The bottom of Fig. 13 shows the time evolution of global attenuation for the algorithms of Table 3 (PR constraint on H P AA ). The I+P adaptation algorithm gives the best results. The transient improvement provided by the various algorithms with respect to the gradient algorithm is slightly less significant than for the algorithms of Table 2. However these algorithms will tolerate high values of the adaption/learning rate.

Broad band disturbance 70-220 Hz

The top of Fig. 14 shows the time evolution of global attenuation for the algorithms given in Table 2. One can draw the same conclusions as for the previous type of disturbance. A significant reduction of the adaptation transient with respect to the gradient is obtained with the various algorithms proposed. The ARIMA2 (2) algorithm (last row of Table 2) provides the best result followed closely by the Conjugate Gradient/Nesterov algorithm. The bottom of Fig. 14 shows the time evolution of global attenuation for the algorithms of Table 3 (PR constraint on H P AA ). Similar results to those obtained for the previous type of disturbance are obtained. The I+P adaptation algorithm gives the best results in terms of adaptation transients.

Conclusions

The paper has tried to give a global picture of the use of dynamic adaptation gain/learning rate for improving the convergence speed of adaptive feedforward noise attenuation schemes using the gradient algorithm with constant adaptation gain/learning rate. The dynamic adaptation gain introduces a frequency weighting which allows to tune it for getting high gain in low frequencies and attenuation in high frequencies. It is this behaviour which explains the benefit of using the dynamic adaptation gain for improving both adaptation transients and steady state performance. While the experimental verifications have been carried out in the context of active noise control, similar results can certainly be obtained in the context of active vibration control.

Appendix A. Average gain property for strictly positive real DAG

Lets us consider the function log C(z -1 )

, where

One seeks to evaluate (43).

Appendix B. Proof of SPR conditions for H DAG

In order to assess the strict real positivity of H DAG (z) on must check the condition

Set z = e jω = cos(ω) + j sin(ω), and the condition (B.1) becomes

f has a finite maximum, and it is located at

. If X max > 1 one must verify f (-1) > 0, moreover one has f (1) > f (-1). If X max < -1 one must verify f (1) > 0, moreover one has f (-1) > f (1). If -1 < X max < 1 one must verify at the same time f (-1) > 0 and f (1) > 0. In any case one must check that min(f (-1), f (1)) > 0. But f (1) > 0 implies that c 1 > -c 2 -1, and f (-1) > 0 implies that c 1 < c 2 + 1. Thus for c 2 < 0 the passivity condition is equivalent to

In this case f is represented by a line, and one must again verify that f (-1) > 0 and f (1) > 0 that leads to the passivity condition -1 < c 1 < 1

In this case f has a finite minimum at

. A sufficient condition for f (X) ≥ 0 ∀X is that f (X) = 0 has a unique solution. In such a situation the discriminant of f denoted ∆ is given by ∆

, and one must have ∆ = 0, which is equivalent to

Thus, one looks for the solutions of (B.3). The discriminant

1 ), and the two solutions of (B.3) are

1 )

On the other hand if -1 ≤ X min ≤ 1 one must have (owing to the expression of X min )

), otherwise this upper bound is given by c 1 < 1 + c 2 , and similarly if c * 1-meets (B.4) the lower bound on c 1 is d and 1 -d ′ 1 > 0. That leads to 1 + c 1 + c 2 < 0. But from the Jury criterion [START_REF] Jury | A stability test for linear discrete time systems in table form[END_REF], a necessary and sufficient condition for 1+c 1 q -1 +c 2 q -2 be a stable polynomial is that at the same time the three following conditions |c 2 | < 1, 1 + c 1 + c 2 > 0, 1 -c 1 + c 2 > 0 are verified. Therefore if β-δ β > 1, one has 1 + c 1 + c 2 < 0 and 1 + c 1 q -1 + c 2 q -2 cannot be stable, and H P AA cannot be positive real. Similarly if β-δ β ≤ -1, one has δ ≥ 2 and since 1 -d ′ 1 > 0 there exists some values of d ′ 1 such that 1 + c 1 + c 2 > 2: for c 1 = 0 this implies c 2 > 1 which is not compatible with the first condition of the Jury criterion, thus in this case 1 + c 1 q -1 + c 2 q -2 cannot be stable and H P AA cannot be positive real. This ends the proof.