Katrin Sieb

Katrin Schuhen

David Sieb

Hubert Wadepohl

Markus Enders
email: markus.enders@uni-hd.de

New Donor Functionalized Cp Ligands: Synthesis and Complexation Behaviour of Quinoxalyl and Benzothiadiazolyl

Keywords: Donor functionalized Cp ligands, Benzothiadiazole, Ferrocene, Chromium, Olefin Polymerization Donor functionalized Cp ligands, Benzothiadiazole, Quinoxaline, Ferrocene, paramagnetic NMR, Olefin Polymerization

Sodium cyclopentadienide reacts as nucleophile with 4,7-dibromo-2,1,3-benzothiadiazole (BTZ) and leads to the new donorfunctionalized ligand Cp BTZ . Related quinoxalyl Cp-systems have been prepared using Pd-catalyzed coupling with zincated Cp-metal complexes. The new ligands comprise two N-donor atoms; one of them is located in a distal position relative to the metal centre so that it cannot coordinate in a chelating manner.. With Cp BTZ ligand derivatives several metal complexes have been synthesized The new chromium(III) complex Cp BTZ CrCl 2 (12) becomes upon activation an active catalyst for the polymerization of ethylene. Relying on DFT-calculations and analysis of spin-density distribution combined with paramagnetic NMR-data a chelating coordination of the Cp BTZ ligand is feasible in 12.

INTRODUCTION

The attachment of neutral donor functions to cyclopentadienes via suitable spacers leads to monoanionic chelating ligands where the steric and electronic properties of cyclopentadienyl ligands are considerably modulated. Therefore a large number of such ligands have been developed 1 and the resulting metal complexes served as active catalysts for the polymerization 2 and oligomerization 3 of olefins, for CH-activation 4 or hydroamination of alkynes 5 . The role of the donor group is the completion of the coordination sphere of the metal in a predefined geometry and to modulate the metal reactivity. Additionally, ligands have been developed where for steric reasons the attached donor group cannot coordinate to the Cp bound metal atom. 6 We now report the synthesis of a new type of donor functionalized Cp ligands where one donor atom is able to coordinate to the Cp bound metal atom and a second donor function remains uncoordinated so that additional metal atoms or Lewis acids can bind. Such interactions could lead to a modulation of the properties of catalysts (e.g. Al-alkyls in olefin polymerization catalysts). The two donor atoms are incorporated into a heteroaromatic ring which enables electronic communication. We have chosen 2,1,3-benzothiadiazole and quinoxaline respectively as pendant groups attached to the Cp rings. This leads to a rigid arrangement of the new ligands with a C 2 spacer between the chelating donor atom and the Cp ring.

RESULTS AND DISCUSSION

Synthesis of Cp-ligands with a pendant 2,1,3-benzothiadiazole

Benzene derivatives usually do not react with cyclopentadienides as nucleophiles except when they are substituted with strongly electron withdrawing groups. 1h, 7 In benzothiadiazoles the homocyclic ring is also activated for the reaction with nucleophiles which can be explained by the contribution of several mesomeric resonance structures (see scheme 1). 8 Therefore we treated 4,7dibromo-2,1,3-benzothiadiazole with lithium indenide (Li-Ind) and after hydrolysis we obtained the donor functionalized indene H1 as a single isomer. The yield of the reaction increases considerably to 90% when two equivalents of lithium indenide are used, as the initially formed product is quickly deprotonated by Li-Ind under formation of the intensely colored anion 1 -. The structure of H1 was established by X-Ray diffraction analysis of a twinned crystal. The four independent molecules found in the asymmetric unit only differ by the relative arrangement of the two ring systems with respect to torsion around the connecting carbon carbon bond (torsion angles -42.6, -43.7, 39.3 and 43.0 °). The distances C1-C2, C1-C6 and C5-C6 are with 1.42 -1.44 Å considerably longer than the other C-C distances in the benzothiadiazole moiety which indicates an important contribution of the mesomeric structure B (scheme 1). 5), C1-C2 1.418( 5), C2-C3 1.365( 5), C3-C4 1.396( 5), C4-C5 1.367( 5), C5-C6 1.437( 5), C6-C1 1.444( 5), C5-C7 1.464( 5), C6-N2-S1 106.2(3), N2-S1-N1 101.5(2), S1-N1-C1 106.3(2), C6-C5-C7-C8 -43.7 [START_REF]Many systems fall into this category. Often the Cp bound metal is already coordinatively saturated. The most prominent examples are donor functionalized ferrocenes. See for example[END_REF].

N S N N S N N S N + 2 R 1 R 2 N S N R 1 R 2 N 1 S N 3 M + M + Br Br Br Br Br Br Br Br -MBr, -HCp R R 1 /R 2 = benzo: 1Li R 1 =R 2 =H: 2Na R 1 = H, R 2 = t-Bu: 3Li H1 N S N t-Bu
The synthesis of H1 shows that lithium indenide is sufficiently reactive for the nucleophilic substitution of bromide in dibromobenzothiadiazole. Compared to lithium indenide, alkaline metal cyclopentadienides are less nucleophilic. Nevertheless the reaction with NaC 5 H 5 leads to the formation of 2Na. After aqueous workup the product could be identified by GC-MS analysis, but it was not possible to isolate the cyclopentadiene H2. The low stability of H2 is in accordance with similar properties of 8quinolylcyclopentadiene. 9 The stability of functionalized cyclopentadienes increases with additional alkyl substituents. Therefore we tried to synthesize the corresponding tetramethyl derivative. However, the unsubstituted carbon atom in lithium tetramethylcyclopentadienide is less nucleophilic compared to the neighboring C-atoms which carry a methyl substituent. Consequently the reaction occurs at a position which is already substituted so that a cyclopentadiene is formed which cannot be deprotonated and hence is no longer useful as ligand (scheme 2). By increasing the steric bulk at the cyclopentadiene substituent a second substitution at this position can be prevented. Therefore lithium t-butylcyclopentadienide was used leading to the isomers H3a and H3b (see scheme 1).

We also tried to attach a second cyclopentadiene to compound 3 by reaction of dibromobenzothiadiazole with an excess of t-butylcyclopentadienide, but this reaction was not successful. However, such a ligand system can be obtained by Negishi-type Pd catalyzed coupling of zincated Cp-metal complexes. Both, mono-zincated and bis-zincated ferrocene react with one equivalent of 4,7dibromobenzothiadiazole in the presence of a Pd(0) catalyst under formation of compound 5. When an excess of the zincated ferrocene is used, the bimetallic compound 6 is obtained (see scheme 3). The ferrocene units in 6 communicate electronically through the benzothiadiazole unit which is expressed by a peak separation of 100 mV in the cyclic voltammogram. We also envisaged synthesizing ligands where a quinoxaline unit is attached to a cyclopentadiene. This could lead to metal complexes with a similar geometry as 8-quinoline metal complexes. The latter compounds have been extensively studied by us and some of them showed interesting catalytic behaviour. 2d, 4 However, the homocyclic part of quinoxaline is not sufficiently activated for nucleophilic substitution so that the Pd-catalyzed reaction with zincated Cp-metal complexes had to be employed. The complexes 7 and 8 have been obtained in good yields and the solid state molecular structures have been determined. Figure 2 and figure 3 show the relative orientation of the nearly planar quinoxaline units to the Cp rings. In 7 the metals are in a trans-arrangement whereas in 8 they are cis-oriented. In both compounds the heterocycles and the Cp rings form angles of approximately 25°. This small angle allows electronic communication between the metal centres which is expressed in the cyclic voltammogram. 7 shows a peak separation of the two oxidation potentials of 80 mV, a value which is somewhat smaller compared to the benzothiadiazole bridged ferrocene 6. 14), C7-C6-C14-C13 24.6 [START_REF] Frisch | [END_REF].

Complexation of ligands 1 -3

The ligands 1 and 3 can be introduced easily into metal complexes, as the protonated ligands are readily available. Two examples with ligand precursor H1 are presented in scheme 5. Deprotonation of H1 with potassium hydride leads to the indenide derivative K1 which was reacted with [ClRh(C 2 H 4 ) 2 ] 2 . The resulting rhodium(I) bis-ethene complex 10 was isolated in 45% yield. As the Rh atom is coordinatively saturated, no N-Rh interaction is possible. Therefore we tried to remove one C 2 H 4 by irradiation with visible light. With the related quinolyl complexes we were able to identify N-donor stabilized Rh(I) mono-ethene complexes which showed activity in catalytic H/D exchange reactions. 4 The irradiation of 10 also leads to ethene elimination, but unfortunately the reaction is not selective and no products could be identified. The zirconium complex 9 was obtained in a clean reaction between H1 and Zr(NMe) 4 . 15 N NMR demonstrates that there is no coordination of the benzothiadiazole to the Zr atom, as the chemical shifts of the two N-atoms are very similar (δ = 333.7, 336.6 ppm) and almost equal to the 15 N shifts of H1 (335.8, 338.6 ppm). This is in contrast to a related indenyl Zr(NMe 2 ) 3 complex where a pendant pyridyl coordinates to the Zr-atom. However, a relatively long Zr-N-distance indicates that this interaction is weak. 11 Thus it is not unexpected that the pendant donor in 9 does not bind to Zr. Deprotonation of the two isomers H3a and H3b with potassium hydride results in the formation of the intensely blue colored anion K3 which shows an absorption maximum in the visible region at 675 nm (ε = 5400). Subsequent reaction with FeCl 2 yields the ferrocene derivative 11 as a blue solid (scheme 6). Due to the unsymmetric 1,3-disubstitution at the Cp-rings, a 1:1 mixture of the C S -isomer and the racemic chiral C 2 -isomers is obtained. The cyclopentadienide Na2 was prepared in situ by the reaction between NaC 5 H 5 and dibromobenzothiadiazole (see scheme 1). Reaction with CrCl 3 (thf) 3 yields the chromium(III) complex 12.

M + R = H Na2 R = t-Bu K3 H 5 H 6 Fe t-Bu BTZ t-Bu BTZ Fe t-Bu BTZ t-Bu BTZ + + FeCl 2 11a 11b
Scheme 6. Synthesis of the ferrocene derivatives 11 and Cr complex 12.

No crystals of 12, suitable for X-ray analysis have been obtained. The mass spectra (EI and HR-FAB) confirm the molecular formula. Complexes with BTZ as ligands are known, 12 but ligand 2 is not ideally suited for a chelating coordination, as the lone-pair at the BTZ substituent is not pointing directly to the metal centre. This results in a constrained geometry of the complex. An other possibility would be a dimer formation with chloro bridges between two Cr centres or a coordination polymer. We calculated monomeric and dimeric molecular structures for 12 with DFTmethods and correlated the calculated Fermi-contact shifts with the experimental paramagnetic shifts observed in the 1 H-NMR spectrum (δ para = δ exp -δ dia ). This procedure works well for chromium(III) complexes with quinolyl-Cp ligands. 13 The theoretical relation between the Fermi-contact shift (δ con ) and the calculated Fermi contact couplings for a S=3/2 system is: 14

δ con = 2.0051•10 5 ρ αβ
The experimental paramagnetic shift δ para is the sum of δ con and δ dip (δ dip is the dipolar chemical shift contribution). For many H-atoms in paramagnetic organometallic compounds the contact shift dominates so that δ para ≈ δ con . A fitting of δ para with ρ αβ with several Cr(III)-Cp-quinolyl complexes leads to the following modified relation: 13 δ con = 1.8484•10 5 ρ αβ + 6.1 ppm Deviations from this simple correlation are significant only for H atoms directly bound to the Cp-Ring and for one position in the quinoline which corresponds to H 6 in complex 12. The DFToptimized structure for monomeric 12 is depicted in Figure 4. d)). 15 Table 1 gives the experimental 1 H-NMR shifts together with calculated chemical shifts using the equation above.

Test of polymerization activity of Cr-complex 12

Chromium(III) complexes with donorfunctionalized Cp ligands lead to highly active catalysts for the polymerization of olefins. 2c-2g Related Cp-Cr complexes with external donors are also active, but less stable. 16 The catalytic activity for the polymerization of ethylene has been tested for 12 and the previously reported 13 14 and the results are compared with polymerization data using zirconocenedichloride as pre-catalyst under identical conditions.

N Me 3 Si

Cr Cl Cl 13

After activation with MAO (1000 equivalents of Al) 12 polymerizes ethylene at room temperature and atmospheric pressure with medium activity. The activity of 13 is much higher and is comparable to the activity to other chromium complexes with 8-quinolyl-Cp ligands or to Cp 2 ZrCl 2 .

CONCLUSION

We have developed the synthesis of benzothiadiazole (BTZ) and quinoxalyl substituted Cp-ligands where one of the two N-donor atoms is connected to the Cp ring by a rigid C 2 -spacer. This ligand geometry is comparable to the geometry in quinolyl substituted Cp ligands (Cp Q ) developed earlier in our group. In addition, a second N-donor is present in the new ligands which is located in a distal position relative to the metal centre. Due to the 5-membered ring in the BTZ ligand the chelating coordination is less favoured compared to Cp Q -ligands. The chromium complex 12 derived from Cp BTZ is active as ethylene polymerization catalyst with moderate activity and produces polyethylene with high molecular weight.

EXPERIMENTAL PART

General Procedures. All experiments were carried out under an atmosphere of dry argon. Solvents were dried by using standard procedures and distilled prior to use. 4,7-dibromo-2,1,3-benzothiadiazole 17 , di-µchlorotetrakis(ethene)dirhodium(I) 18 , 1,1'-dilithioferrocene 19 were prepared by literature procedures. All other reagents were used as purchased. NMRspectra were recorded on a Bruker DRX 200, Bruker Avance II 400 and Bruker Avance III 600 spectrometer ( 1 H: 200 MHz, 400 MHz and 600 MHz respectively; 13 C: 50 MHz, 100 MHz and 150 MHz respectively). The spectra were referenced using residual proton peaks or the 13 C solvent peaks. 15 N NMR spectra were recorded on a Bruker Avance III 600 spectrometer equipped with a cryogenically cooled probe (QNPcryoprobe TM ). Infrared spectra were obtained on a Varian 3100 FT-IR spectrometer (Excalibur Series). UV/Vis spectra have been measured in 1 cm cells with a Varian Cary 5000 spectrometer. Mass spectra and elemental analysis were performed by the analytical services at the Organisch-Chemisches Institut, University of Heidelberg. Complexes 10and 11 did not give satisfactory elemental analysis data. Cyclic voltammetry measurements were done with a PAR 263A potentiostat, a METROHM glassycarbon-disk working electrode, a saturated calomel reference electrode, and dichloromethane solvent with 0.1 mol L -1 (n-Bu4)NPF6 as supporting electrolyte; the scan rate was 100 mV s -1 . Molecular weight distributions of polyethylene were measured with a Polymer Lab HT-GPC (PL-GPC 220) at 150°C in 1,2,4-trichlorobenzene. Molecular masses are calibrated to PE. 4-Bromo-2,1,3-benzothiadiazol-7-yl-1-indene (H1) A solution of 711 mg (6.12 mmol) indene in 50 ml thf was cooled to -50°C and 3.8 ml (6.12 mmol) n-BuLi (1.6 M in hexane) were added dropwise and stirred at this temperature for 1 hour. At room temp. the solution was added via cannula to a solution of 1.09 g (3.70 mmol) of 4,7-dibromo-2,1,3-benzothiadiazole in 100 ml thf. After 24 hours at room temp. a deep green reaction mixture was obtained. After addition of 100 ml of aqueous HCl (5%) the mixture was treated with NH4OH until a pH of 12 was obtained. The layers were separated and the aqueous one was extracted twice with diethylether (2 x 100 ml). The organic layers were collected, dried over MgSO4, filtered and the volatile components were removed to obtain a yellow powder. After column chromatography (Al2O3) with toluene as eluent and removing of the solvent a yellow powder was obtained. Yield: 1.10 g (3.34 mmol = 90%). 1 H NMR (CDCl3, 200 MHz) δ = 3.65 (d, 2H, 3 J (H, H) = 2.0 Hz, CH2), 7.10 (t, 1H, 3 J(H, H) = 2.0 Hz, H 12 ), 7. 26 -7.32 (m, 2HInd), 7.43 (d, 1H, 3 J(H, H) = 6.3 Hz, HInd), 7.57 (d, 1H, 3 J(H, H) = 6.3 Hz, HInd), 7.69 (d, 1H, 3 J(H, H) = 7.5 Hz, H 6 ), 7.91 (d, 1H, 3 J(H, H) = 7.5 Hz, H 5 ) ppm. 13 C NMR (CDCl3, 50 MHz) δ = 38.9 (CH2), 113.2 (Cq), 120.5 (CH 12 ), 124.3, 125.3, 126.2, 128.1 (4 x CHInd), 128.9 (Cq), 132.1, 136.1 (CH 5 and CH 6 ), 139. 9, 143.3, 144.2, 153.5, 153.8 (5 x Cq) ppm. 15

4-Bromo-2,1,3-benzothiadiazol-7-yl-tert-butylcyclopentadiene (H3)

To a solution of 1.64 g (15.4 mmol) of 6,6-dimethylfulvene in 30 ml of diethylether was added 9.65 ml (15.4 mmol) of methyllithium (1.6 M solution in diethylether) at 0°C. After 5 min the solvent was evaporated and the resulting white solid was washed with diethylether, dried and dissolved in 50 ml of thf. This solution was added to a solution of 2.70 g (9.3 mmol) of 4,7-dibromo-2,1,3-benzothiadiazole in 50 ml of thf within 10 min. The reaction mixture was heated to 40°C for 21 h to give an intensely blue colored solution which was then added to a saturated NH4Cl solution in water. The aqueous phase was extracted twice with n-hexane and the combined organic phases were dried over MgSO4, filtered and the solvent evaporated. The crude product was purified by chromatography (Al2O3, eluent: toluene) to give 925 mg (30%) of the two isomers of H3 in a 1:1 ratio. Isomer H3a: 1 H NMR (CDCl3, 600 MHz) δ = 1.25 (s, 9H, CH3), 3.56 (dd, 2H, 4 J (H 9 ,H 11 ) = 1.3 Hz, 4 J (H 9 ,H 12 ) = 1.3 Hz, H 9 ), 6.34 (m, H 11 ), 7.39 (d, 1H, 3 J(H 6 ,H 5 ) = 7.8 Hz, H 6 ), 7.73 (d, 1H, 3 J(H 5 ,H 6 ) = 6.9 Hz, H 5 ), 7.77 (m, H 12 ) ppm. 13 C NMR (CDCl3, 100 MHz) δ = 30.8 (CH3), 32.3 (Cq), 40.9 (CH2), 110.4 (Cq), 124.2 (CH 6 ), 125.7 (CH 1 ), 132.3 (CH 5 ), 132.4 (Cq), 135.8 (CH 12 ), 139. 1, 152.5, 154.0, 161.5 (4 x Cq) ppm. Isomer H3b: 1 H NMR (CDCl3, 600 MHz) δ = 1.25 (s, 9H, CH3), 3.54 (dd, 2H, 3 J(H 12 ,H 11 ) = 1.5 Hz, 4 J(H 12 ,H 9 ) = 1.5 Hz, H 12 ), 6.12 (td, 1H, 3 J(H 11 ,H 12 ) = 1.5 Hz, 4 J(H 11 ,H 9 ) = 1.6 Hz, H 11 ), 7.44 (d, 1H, 3 J(H 6 ,H 5 ) = 7.8 Hz, H 6 ), 7.75 (d, 1H, 3 J(H 5 ,H 6 ) = 6.9 Hz, H 5 ), 7.99 (dt, 1H, 3 J(H 9 ,H 11 ) = 1.6 Hz, 4 J(H9,H12) = 1.5 Hz, H9) ppm. 13 C-NMR (CDCl3, 100 MHz) δ = 29.7 (CH3), 33.6 (Cq), 41.4 (CH2), 111.0 (Cq), 123.9 (CH 11 ), 124.8 (CH 10 ), 132.3 (CH 5 ), 132.4 (Cq), 136.3 (CH 9 ), 141.9, 152.5, 154.0, 158.8 (4 x Cq) ppm. MS (EI) m/z (%): 334 (53) [M + ], 321 (100) [M + -CH3], 278 (57) [M + -t Bu], 240 (36) [M + -CH3 -Br]. UV/Vis (thf), λmax [nm] (ε [l mol -1 cm -1 ]): 306 (9740), 438 (4550). Deprotonation with KH leads to the intensely colored anion 3K: UV/Vis (thf), λmax [nm] (ε [l mol -1 cm -1 ]): 339 (9700), 675 (5400).

4-Bromo-2,1,3-benzothiadiazol-7-yl-tetramethylcyclopentadiene (4)

A solution of lithium tetramethylcyclopentadienide in 50 ml thf was prepared from 830 mg (6.8 mmol) of 2,3,4,5-tetramethylcyclopentadiene and 4.25 ml (6.8 mmol) of a n-BuLi solution (1.6 M in hexane). This solution was added to a solution of 1.0 g (3.4 mmol) 4,7-dibromo-2,1,3benzothiadiazole in 100 ml thf at -55°C. The reaction mixture turned blue and after stirring at room temp. for 18 h the solution was violet. Diluted hydrochloric acid was added and the mixture was alkalized with NH4OH until a pH of 12 was obtained. The mixture was extracted with diethylether and the combined organic layers were dried with MgSO4. The product mixture was purified by chromatography and the main component was identified as compound 4. Yield 450 mg (1.34 mmol), 39%. 1 H NMR (CDCl3, 400 MHz) δ = 1.68 (s, 3H, CH3), 1.78 (s, 3H, CH3), 1.85 (s, 6H, 2 x CH3), 6.24 (s, 1H, H 12 ), 7.03 (d, 1H, 3 J (H,H) = 7.6, H 6 ), 7.63 (d, 1H, 3 J (H,H) = 7.6, H 5 ) ppm. 13 C-NMR (CDCl3, 100 MHz): δ = 11.0, 11.3, 14.0, 20.5 (4 x CH3), 58.8 (Cq 8 ), 111.8 (Cq), 125.8 (CH 6 ), 132.2 (CH 7 ), 135.7 (Cq), 136.1 (CH 12 ), 136.9, 141.0, 143.8 (3 x Cq-Cp) C15H15BrN2S (335.26) calcd. C 54.74,H 4.51,N 8.36;found C 53.89,H 4.52,N 8.28.

4-Bromo-7-ferrocenyl-2,1,3-benzothiadiazole (5)

To a solution of 1.50 g (8.1 mmol) of ferrocene in 30 ml thf at -40°C was added 5.5 ml (6.7 mmol) of a 1.26 M solution of t-BuLi in pentane. After stirring for 2 h at room temp. a solution of 884 mg (6.5 mmol) of ZnCl2 was added and the resulting solution of ferrocenylzinc chloride was stirred for an additional hour. A solution of Pd(PPh3)2 (prepared from 235 mg (0.35 mmol) of (PPh3)2PdCl2 and 0.70 mmol of DibAl-H) was added to the red solution of ferrocenylzinc chloride and subsequent 1.88 g (6.4 mmol) of 4,7-dibromobenzo-2,1,3-thiadiazole in 50 ml thf was added. After stirring at room temp. for 5 days a deep violett solution was obtained which was treated with a solution of 5 g (0.2 mol) NaOH in 100 ml H2O. The layers were separated and the aqueous layer was extracted twice with CH2Cl2. The organic layers were dried with MgSO4, filtered and the volatile components were removed to obtain a violet powder with ferrocene impurity. Purification by column chromatography (Al2O3 (neutral); eluent: toluene) led after removing of the solvent to violet crystals. Yield: 1.02 g (2.56 mmol), 40%. 1 H NMR (CDCl3, 400 MHz) δ = 3.99 (s, 5H, Cp-H), 4.48 (pt, 2H, Cp-CH), 5.22 (pt, 2H, Cp-CH), 7.52 (d, 1H, 3 J(H,H) = 7.7 Hz, H 6 ), 7.70 (d, 1H, 3 J(H,H) = 7.7 Hz, H 5 ) ppm. 13 C NMR (CDCl3, 150 MHz) δ = 68.9 (CH(C5H5), 70.0, 70. 1 (2 x CH(C5H4), 68.3 (Cq), 113.7 (Cq), 132.3, 132.6 (CH 5 and CH 6 ), 133.9, 152.8, 152.9 (3 x Cq) ppm. MS (EI) m/z (%): 398 (M + , 100), 333 (M + -Cp, 25), 277 (M + -FeCp, 18). MS (HR-FAB) calcd. 397.9176; found 397.9180.

4,7-(Bis-ferrocenyl)-2,1,3-benzothiadiazole (6)

A solution of 1,1'-bis-chlorozinc-ferrocene in 100 ml thf was prepared from 3.11 g (11.7 mmol) 1,1'-dilithioferrocene x 2/3 TMEDA and 3.18 g (23.3 mmol) ZnCl2. A solution of (PPh3)2Pd (prepared from 170 mg (0.23 mmol) of (PPh3)2PdCl2 and 0.48 mmol of DibAl-H) and 2.0 g (6.8 mmol) of 4,7dibromo-2,1,3-benzothiadiazole in 20 ml thf were added slowly to the reaction mixture leading to a deep red solution. After 3 days stirring at room temp. 10 g (0.23 mol) NaOH in 100 ml H2O were added. After extraction with CH2Cl2 (2 x 100 ml) the organic layers were dried with MgSO4 and the solvent was removed. Column chromatography (Al2O3 (neutral), eluent: toluene) led to a violet powder. Yield: 0.81 g (1.61 mmol), 24 %. 1 H NMR (C6D6, 200 MHz) δ = 3.95 (s, 10H, Cp-H), 4.33 (pt, 4H, Cp-H), 5.33 (pt, 4H, Cp-H), 7.41 (s, 2H, CH) ppm. 13 C NMR (C6D6, 50 MHz) δ = 68.9, 69.9 (2 x CH(C5H4)), 70.1 (CH(C5H5)), 82.2 (Cq-Cp), 125.5 (CH), 131.2, 154.6 (2 x Cq) ppm. -MS (FAB + ) m/z (%) 504 (100) [M + ]. C26H20N2SFe2 (504.21) calcd. C 61.93,H 4.00,N 5.56;found C 61.87,H 4.47,N 5.18. CV: E 0/+1 = 0.455 V and = 0.555 V vs. SCE.

5,8-(Bis-ferrocenyl)-quinoxaline (7)

A solution of 1,1'-(bis-chlorozinc)ferrocene in 100 ml thf was prepared from 2.01 g (7.55 mmol) 1,1'-dilithioferrocene x 2/3 TMEDA and 2.04 g (15.1 mmol) ZnCl2. A solution of (PPh3)2Pd (prepared from 265 mg (0.38 mmol) of (PPh3)2PdCl2 and 0.76 mmol of DibAl-H) and 0.55 g (1.9 mmol) of 4,7-dibromo-2,1,3-benzothiadiazole in 20 ml thf were added slowly to the reaction mixture leading to a red-violet solution. After 5 days stirring at room temp. 10 g (0.23 mol) NaOH in 100 ml H2O were added. After extraction with CH2Cl2 (2 x 100 ml) the organic layers were dried with MgSO4 and the solvent was removed. Column chromatography (Al2O3 (neutral), eluent: toluene) led to red-violet crystals. Yield: 0.33 g (0.66 mmol), 35%. 1 H NMR (CDCl3, 200 MHz) δ = 4.09 (s, 10H, Cp-H), 4.46 (pt, 4H, Cp-CH), 5.18 (pt, 4H, Cp-CH), 7.99 (s, 2H, H 6 and H 7 ), 8.86 (s, 2H, H 2 and H 3 ) ppm. 13 A solution of (ClZnC5H4)Re(CO)3 was prepared by addition of 5.4 ml (4.0 mmol) of n-BuLi (1.6 M in hexane) to a solution of 1.21 g (3.6 mmol) of C5H5Re(CO)3 in 70 ml of thf at -78°C. After 1h at this temperature a solution of ZnCl2 in 70 ml of thf was added. 1h later, a solution of (PPh3)2Pd (prepared from 140 mg (0.2 mmol) of (PPh3)2PdCl2 and 0.4 ml (0.4 mmol) of DibAl-H (1M in hexane)) and subsequently a solution of 1.04 g (3.6 mmol) of 5,8-dibromoquinoxaline in 20 ml of thf was added. The reaction mixture was allowed to warm to room temp. and after stirring for 5 days 100 ml of NaOH solution was added. The product was extracted with CH2Cl2 and purified by chromatography (Al2O3/toluene) to give 1.86 g (2.34 mmol, 65%) of yellow solid 8. 1 H NMR (CDCl3, 400 MHz) δ = 5.46 (pt, 4H, CHCp), 6.32 (pt, 4H, CHCp), 7.81 (s, 2H), 8.90 (s, 2H) ppm. 13 C-NMR: (CDCl3, 100 MHz): δ = 84. 4,88.3 (CHCp),128.5 (CH),131.5,140.5 (2 x Cq),143.8 (CH),194.2 (CO) 9) In a teflon valve NMR-tube 15 mg (46 µmol) of H1 was dissolved in C6D6 and 12 mg (46 µmol) of tetrakis(dimethylamido)zirconium(IV) was added. The solution turned into deep red. NMR spectra were recorded and showed compound 9 together with one equivalent of HNMe2. Yield: > 90% (not isolated). 1 H NMR (C6D6, 600 MHz) δ = 2.69 (s, 18H, 3 x N(CH3)2), 6.56 (dd, 1H, 3 J(H 11 ,H 12 ) = 3.6 Hz, 4 J(H 11 ,HInd) = 0.6 Hz, CH 11 ), 7.09 (ddd, 1H, 3 J = 8.3 Hz, 3 J = 6.5 Hz, 4 J = 0.9 Hz, HInd), 7.16 (ddd, 3 J = 8.3 Hz, 3 J = 6.5 Hz, 4 J = 0.9 Hz 1H, H3´), 7.56 (d, 1H, 3 J = 7.7 Hz, CH 5 or CH 6 ), 7.59 (d, 1H, 3 J = 7.7 Hz, CH 5 or CH 6 ), 7.67 (ddd, 1H, 3 J = 8.3 Hz, 4 J = 0.9 Hz, 5 J = 0.9 Hz, HInd), 7.71 (d, 1H, 3 J(H 11 ,H 12 ) = 3.6 Hz, CH 12 ), 7.99 (dddd, 1H, 3 J = 8.3 Hz, 4 J = 0.6 Hz, 4 J = 0.9, 5 J = 0.9 Hz, HInd) ppm. 13 C NMR (C6D6, 150 MHz) δ = 43.8 (N(CH3)2), 99.1 (CH 11 ), 108.3, 110.5 (2 x Cq), 119.8 (CH 12 ), 121.8, 123.0 (2 x CHInd), 123.6 (Cq), 123.7, 123.8 (2 x CHInd), 126.2 (CH 5 or CH 6 ), 127.4, 129.9 (Cq), 132.3 (CH 5 or CH 6 ), 153.6, 154.1 (2 x Cq) ppm. 15 N NMR (C6D6, 61 MHz): δ = 333.7, 336.6 ppm. η 5 -(4-Bromo-2,1,3-benzothiadiazol-7-yl-indenyl)-bis(η 2 -ethen)]rhodium (10) In a schlenk flask 99 mg (0.30 mmol) of H1 was dissolved in 20 ml thf and 12 mg (0.30 mmol) KH was added. After 12 hours stirring at room temp. the solution turned green. A solution of 58 mg (0.15 mmol) [ClRh(C2H4)2]2 in 20 ml toluene was added dropwise. The reaction mixture was stirred for 24 hours and became red. After removing of the solvent in vacuum a red powder was obtained. After extraction with toluene and filtration of non-soluble residues, the product was purified by column chromatography (Al2O3 (neutral), eluent: toluene). The product was obtained as a red powder. Yield: 66 mg (0.13 mmol), 45%. 1 H NMR (CDCl3, 400 MHz) δ = 1.90, 2.22 (br, 8H, CH2), 5.18 (d, 1H, 3 J(H,H) = 2.4 Hz, HCp); 7. 24 -7.27 (m, 3H, 2 x HInd and HCp), 7.46 (d, 1H, 3 J(H,H) = 7.3 Hz, HInd), 7.52 (d, 1H, 3 J(H,H) = 7.9 Hz, H 6 ), 7.63 (d, 1H, 3 J(H,H) = 7.3, HInd), 7.88 (d, 1H, 3 J(H,H) = 7.9 Hz, H 5 ) ppm. 13 C NMR (CDCl3, 100 MHz) δ = 47.7 (d, 1 J(Rh,C) = 13.7 Hz, C2H4), 78.2 (d, 1 J(Rh,C) = 4.9 Hz, CHCp), 95.1 (d, 1 J(Rh,C) = 5.2 Hz, CHCp), 110.9, 111.6 (2 x Cq), 112.5 (d, 1 J(Rh,C) = 2.6 Hz, Cq), 118.3,120.6,123.9,124.5 (4 x CHInd)

Bis-[η 5 -(4-bromo-2,1,3-benzothiadiazol-7-yl)-cyclopentadienyl]iron (11)

To a solution of 183 mg (0.54 mmol) 4-bromo-2,1,3-benzothiadiazol-7-yl-tbutylcyclopentadiene in 10 ml thf was added 22 mg (0.55 mmol) of potassium hydride. After 10 h stirring at room temperature this solution is dropped into a suspension of 40 mg (0.32 mmol) FeCl2 in 10 ml of thf. After 2 d the solvent is evaporated and the crude product is purified by chromatography (Al2O3, toluene:hexane = 2:1). Yield 219 mg (0.3 mmol), 56%. 1 H NMR (CDCl3, 600 MHz) δ =1.01, 1.24 (2 x s, 2 x 9H, CH3), 4.59,4.75,4.67,5.01,6.56,7.44,. MS (EI): m/z (%): 724 [START_REF] Beurskens | Crystallographic Computing[END_REF] [M+], 644 (8) [M+ -Br], 566 (13) [M+-2 x Br]. UV/Vis (thf), λmax [nm] (ε) [l mol -1 cm -1 ]): 290 (76000), 388 (36000), 562 (11000).

[η 5 -(4-Bromo-2,1,3-benzothiadiazol-7-yl)-cyclopentadienyl]chromium(III)dichloride (12)

A solution of 555 mg (6.3 mmol) sodium cyclopentadienide in 60 ml thf was added dropwise to 929 mg (3.16 mmol) 4,7-dibromo-2,1,3benzothiadiazole in 70 ml thf. The reaction mixture was stirred at room temp. for 24 h leading to a deep blue color. This solution was added slowly via a cannula to 1.18 g (3.16 mmol) CrCl3(thf)3 in 80 ml of thf. After 2 days at room temp. the reaction mixture turned to intensely green. Removal of the solvent led to a green powder which was washed with hexane and the product was extracted with toluene. After evaporation of the solvent 490 mg (1.22 mmol, 39%) of the solid green complex 11 were isolated. 1 conventional direct methods 20 , 21 (H1 and 8) or heavy atom method combined with structure expansion by direct methods 22 (7). Refinement: full-matrix least squares methods based on F 2 ; all non-hydrogen atoms anisotropic. Hydrogen atoms were input at calculated positions and refined riding (H1 and 8) or fully refined (7). 23 , 21 The twin refinements were against all overlapping and non-overlapping reflections. Crystallographic data (excluding structure factors) for the structures reported in this paper have been deposited with the Cambridge Crystallographic Data Center: CCDC-713319 (H1), -713317 (7) and -713318 (8). These data can be obtained free of charge at http://www.ccdc.cam.ac.uk/cgi-bin/catreq.cgi? [or from CCDC, 12 Union Road, Cambridge CB2 1EZ, UK; Phone: +44 1223 336408, Fax: +44 1223 336033; E-mail: deposit@ccdc.cam.

Scheme 1 :

 1 Scheme 1: Mesomeric resonance structures of benzothiadiazole (A -C) and formation of cyclopentadienes functionalized with 2,1,3benzothiadiazole. (Numbering of atoms for NMR assignment in exp. part).

Figure 1 .

 1 Figure 1. Solid-state molecular structure of one of four independent molecules of H1. Selected bond lengths [Å] and angles [°]: C1-N1 1.347(5), N1-S1 1.617(3), S1-N2 1.611(3), N2-C6 1.357(5), C1-C2 1.418(5), C2-C3 1.365(5), C3-C4 1.396(5), C4-C5 1.367(5), C5-C6 1.437(5), C6-C1 1.444(5), C5-C7 1.464(5), C6-N2-S1 106.2(3), N2-S1-N1 101.5(2), S1-N1-C1 106.3(2), C6-C5-C7-C8 -43.7[START_REF]Many systems fall into this category. Often the Cp bound metal is already coordinatively saturated. The most prominent examples are donor functionalized ferrocenes. See for example[END_REF].

  Scheme 2. Formation of the "wrong" isomer 4.

  10

  Scheme 3. Formation of mono-and bis-ferrocenyl benzothiadiazoles.

8 2 4 .

 84 Scheme 4. Formation of quinoxalyl bridged complexes 7 and 8.

Figure 2 .

 2 Figure 2. Solid state molecular structures of 7. Selected bond lengths [Å] and angles [°]: C5-C6 1.478(2), C6-C7 1.429(3), C7-N1 1.357(2), N1-C8 1.318(3), C8-C9 1.406(3), C9-N2 1.321(3), N2-C10 1.368(2), C10-C11

  Scheme 5. Complexation of ligand 1.

Figure 4 .

 4 Figure 4. DFT-optimized structure of 12 (UB3LYP/6-311g(d)).15

  , 154.0, 154.4 (2 x Cq) ppm. MS (EI) m/z (%): 334 (80) [M + ], 319 (50) [M + -CH3], 240 (100) [M + -Br-CH3], 225(85) [M + -C2H6 -Br].

  C NMR (CDCl3, 50 MHz) δ = 69.1 (CH(C5H4)), 69.4 (CH(C5H5)), 70.5 (CH(C5H4)), 82.6 (Cq-Cp), 128.7 (C 6 and C 7 ), 136.5, 141.4 (2 x Cq), 142.9 (C 2 and C 3 ) ppm. MS (FAB + ) m/z (%) 498 (100) [M + ], 433 (25) [M + -Cp], 368 (18) [M + -2 Cp]. MS (HR-FAB) calcd. 498.0481; found 498.0470. C28H22Fe2N2 (498.05): calcd.C 67.51, H 4.45, N 5.62; found C 67.96, H 4.63, N 5.40. CV: E 0/+1 = 0.415 V and 0.495 V vs. SCE. η 5 -(5,8-Bis-cyclopentadienylchinoxalyl)]rhenium(I)(8)

  ppm. MS (EI) m/z (%): 796 (10) [M + ], 768 (40) [M + -CO], 740 (40) [M + -2 CO], 712 (38) [M + -3 CO], 684 (25) [M + -4 CO], 656 (20) [M + -5 CO], 628 (20) [M + -6 CO]. FT-IR (toluene) ν(tilde) (cm -1 ) = 2021 (s, νC=O), 1928 (s, νC=O). C24H12N2O6Re2 (796.78): calcd. C 36.18, H 1.52, N 3.52; found C 36.12

  H NMR (C6D6, 400 MHz) δ = 159 (s, ν½ ~ 4000 Hz, CHCp), 19 (s, ν½ ~ 500 Hz, CH 5 ), 6.8 (s, ν½ ~ 18 Hz, CH 6 ) ppm. MS (EI) m/z (%): 364 (25) [M + -Cl], 329 (100) [M+ -2 Cl]. MS (HR-FAB) C11H6N2 79 Br 35 ClSCr calcd. 363.8529; found 363.8555. Crystal structure determinations Crystal data. H1: C15H9BrN2S, monoclinic, space group P21, a = 19.1008(14), b = 7.4055(6), c = 19.1311(14) Å, β = 109.4620(10) °, V = 2551.5(3) Å 3 , Z = 8, µ = 3.370 mm -1 , F000 = 1312. T = 100(2) K, θ range 1.1 to 30.1 °. Index ranges h, k, l (indep. set): -26 ... 26, -10 ... 10, -26 ... 26. Reflections measd.: 59944, indep.: 14928 [R int = 0.0448], obsvd. [I > 2σ(I)]: 13768. Final R indices [F o > 4σ(F o )]: R(F) = 0.0313, wR(F 2 ) = 0.0620, GooF = 1.018. Pseudo-merohedral and racemic twin. 7: C28H22Fe2N2, monoclinic, space group P21/c, a =14.2530(16) , b = 7.4496(8), c = 22.8163(19) Å, β = 123.218(5) °, V = 2026.7(4) Å 3 , Z = 4, µ = 1.452 mm -1 , F000 = 1024. T = 100(2) K, θ range 1.7 to 32.0°. Index ranges h, k, l (indep. set): -21 ... 20, -11 ... 0, -34 ... 19. Reflections measd.: 34112, indep.: 6945 [R int = 0.0645], obsvd. [I > 2σ(I)]: 5121. Final R indices [F o > 4σ(F o )]: R(F) = 0.0410, wR(F 2 ) = 0.1081, GooF = 1.058. 8: C24H12N2O6Re2, triclinic, space group P-1, a = 8.0187(16), b = 9.5933(19), c = 15.026(3) Å, α = 71.948(4), β = 123.218(5), γ = 82.969(4) °, V = 1072.8(4) Å 3 , Z = 2, µ = 2.467 mm -1 , F000 = 736. T = 150(2) K, θ range 2.2 to 32.2°. Index ranges h, k, l (indep. set): -11 ... 11, -13 ... 14, 0 ... 22. Reflections measd.: 28868, indep.: 7906 [R int = 0.0839], obsvd. [I > 2σ(I)]: 5934. Final R indices [F o > 4σ(F o )]: R(F) = 0.0703, wR(F 2 ) = 0.1096, GooF = 1.002. Pseudo-merohedral twin. Intensity data collection: Bruker AXS Smart 1000 CCD diffractometer, Mo-Kα radiation, graphite monochromator, λ = 0.71073 Å. Structure solution:

  Figure 1 241x137mm (600 x 600 DPI)

Table 1 :

 1 Comparison of experimental 1 H-NMR shifts with calculated shifts.

	H-atom	δ(exp)	δ(para)	ραβ [a.u. x 10 5 ]	δcon
	H-Cp	159	153	72 / 98	139,2 / 187,2
	H 5	19	11	12	28,3
	H 6	6.8	-3	-18	-27,2

Table 2 :

 2 Polymerization results a)

	Pre-catalyst	activity b) [g mmol -1 h -1 ]	MW	PDI
	13	2550	247.000	2,7
	Cp2ZrCl2	2330	607.000	2,4
	12	330	824.000	3.6

a) Polymerization conditions: 100 ml toluene, 25°C, 3 µmol Kat., 1000 eq. MAO (calcd. for ratio Al:Cr(Zr)), ethene pressure 1 bar, 20 min. b) Under the conditions used, activities greater 1500 become more and more levelled, as diffusion limitation leads to lower olefin concentrations.

Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (SFB 623).

Wiley-VCH ZAAC


Olivier Bouissou
email: olivier.bouissou@cea.fr

Yassamine Seladji
email: yassamine.seladji@cea.fr

Alexandre Chapoutot
email: alexandre.chapoutot@lip6.fr

Abstract Fixpoint Computations with Numerical Acceleration Methods *

Keywords: Abstract numerical domains, acceleration of convergence, widening operator

come   L'archive ouverte pluridisciplinaire

Introduction

In the field of static analysis of embedded, numerical programs, abstract interpretation [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF][START_REF] Cousot | Comparing the Galois connection and widening/narrowing approaches to abstract interpretation[END_REF] is widely used to compute over-approximations of the set of behaviors of a program. This set is usually defined as the least fixpoint of a monotone map on an abstract domain given by the (abstract) semantics of the program. Using Tarski's theorem [START_REF] Tarski | A lattice-theoretical fixpoint theorem and its applications[END_REF], this fixpoint is computed as the limit of the iterates of the abstract function starting from the least element. These iterates build a sequence of abstract elements that (order theoretically) converge towards the least fixpoint. This sequence converging often slowly (or even after infinitely many steps), the theory of abstract interpretation introduces the concept of widening [START_REF] Cousot | Comparing the Galois connection and widening/narrowing approaches to abstract interpretation[END_REF].

A widening operator is a two-arguments function ∇ which tries to predict the limit of the iterates based on the relative position of two consecutive iterates. For example, the standard widening operator on the interval abstract domain consists in comparing the limits of the intervals and setting the unstable ones to ∞ (or -∞). A widening operator often makes large over-approximation because it must make the sequence of iterates converge in a finite time. Over-approximation may be reduced afterward using a narrowing operator but the precision of the final approximation still strongly depends on the precision of the ∇. Various techniques have been proposed to improve it. Delayed widening makes use of ∇ after n iteration steps only (where n is a user-defined integer), thus letting the first loop iterates execute before trying to predict the limit. Another approach is to use a widening with thresholds [START_REF] Blanchet | Design and implementation of a special-purpose static program analyzer for safety-critical real-time embedded software[END_REF]: the upper bound of the interval (for example) is not directly set to ∞, but is successively increased using a set of thresholds that are candidates for the value of the fixpoint upper bound. In practice, these techniques are necessary to obtain precise fixpoint approximations for industrial sized embedded programs. However, they suffer from their lack of automatization: thresholds must be chosen a priori and are defined by the user (to the best of our knowledge, no methods exist to automatically find the best thresholds). The delay parameter n is also to be defined a priori. This makes the use of a static analyzer difficult as these (non trivial) parameters are often hard to find.

In this article, we present some ongoing work which shows that it is possible to use sequence transformation techniques in order to automatically and efficiently derive approximation of the limit of Kleene iterates. This approximation may not be safe (i.e. may not contain the actual limit), but we show how to use it in the theory of abstract interpretation. Sequence transformation techniques (also known as convergence acceleration methods) are widely studied in the field of numerical analysis [START_REF] Brezinski | Extrapolation Methods-Theory and Practice[END_REF]. They transform a converging sequence (x n ) n∈AE of real numbers into a new sequence (y n ) n∈AE which converges faster to the same limit (see Section 3.2). In some cases (depending on the method), the acceleration is such that (y n ) n∈AE is ultimately constant. Some recent work [START_REF] Brezinski | A review of vector convergence acceleration methods, with applications to linear algebra problems[END_REF] applied these techniques in the case of sequences of vectors of real numbers: vector sequence transformations introduce relations between elements of the vectors and perform better than scalar ones. Our main contribution is to show that we can use these methods in order to improve the fixpoint computation in static analysis: we define dynamic thresholds for widening that are very close to the actual fixpoint. This increased precision is obtained because the sequence transformations use all iterates and quantitative information (i.e. relative to the distance between elements) to predict the limit. They thus have access to more information than the widening operator and can make better prediction. In this work, we focus on the interval domain, but we believe that this work may be applied for any abstract domain, especially the ones with a pre-defined shape (octagons [START_REF] Miné | Weakly Relational Numerical Abstract Domains[END_REF], templates [START_REF] Sankaranarayanan | Scalable analysis of linear systems using mathematical programming[END_REF], etc.).

This article is organized as follows. In Section 2, we explain on a simple example how acceleration methods may be used to speed-up the fixpoint computation. In Section 3, we recall the theoretical basis of this work and present our main theoretical contribution. Section 4 presents some early experiments on various floating-point programs that show the interest of our approach, while Sections 5 and 6 discuss related works and perspectives.

Notations. In the rest of this article, (x n ) will denote a sequence of real numbers (i.e. (x n ) ∈ Ê AE ), while (x n ) denotes a sequence of vector of real numbers (i.e. (x n ) ∈ Ê p AE for some p ∈ AE). The symbol X n will be used to design abstract iterates, i.e. X n ∈ A for some abstract lattice A.

An introductive example

In this section, we explain, using a simple example, how sequence acceleration techniques can be used in the context of static analysis. In short, our method works as follows: let (X n ) be a sequence of intervals computed by the Kleene iteration and that is chosen to be widened (see [START_REF] Bourdoncle | Efficient chaotic iteration strategies with widenings[END_REF] for details on how to chose the widening points). We extract from (X n ) a vector sequence (x n ): at stage k, x k is a vector that contains the infimum and supremum of each variable of the program. As Kleene iterates converge towards the least fixpoint of the abstract transfer function, the sequence (x n ) converges towards a limit x which is the vector containing the infimum and supremum of this fixpoint. We then compute an accelerated sequence (y n ) that converges towards x faster than (x n ). Once this sequence has reached its limit (or is sufficiently close to it), we use x as a threshold for a widening on (x n ) and thus obtain, in a few steps, the least fixpoint. In the rest of this section, we detail these steps.

The program. We consider a linear program which iterates the function F (X) = A • X + B • U where A, B and U are constant matrices and X is the vector of variables (see Figure 1). Initially, we have

x1 ∈ [1, 2], x2 ∈ [1, 4], x3 ∈ [1, 20], u1 ∈ [1, 6], u2 ∈ [1, 4] and u3 ∈ [1, 2]
. Using an existing analyzer working on the interval abstract domain, we showed that this program converges in 55 iterations (without widening) and obtained the invariant [-5.1975, 8.8733] for x1 at line 2.

Extracting the sequences. From this program, we can define a vector sequence of size 6,

x n = x 1 n , x 1 n , x 2 n , x 2 n , x 3 n , x 3
n , which represent the evolution of the suprema and infima of the variables x1, x2 and x3 at line 2. For example, the sequence (x 1 n ) is recursively defined by: Note that we are not interested in the formal definition of these sequences (as given by Equation ( 1)), but only in their numerical values that are easily extracted from Kleene iterates. Each sequence (x i n ) (resp. (x i n )) is increasing (resp. decreasing) and the sequence (x n ) converge towards a vector x containing the infimum and supremum of the fixpoint (see Figure 2, dotted lines).

x 1 n+1 = max x 1 n , -0.4375 * x 1 n + 0.0625 * x 2 n + 0.2652 * x 3 n + 0.1 * u1 . (1
Accelerating the sequences. We then used the vector ε-algorithm [START_REF] Brezinski | A review of vector convergence acceleration methods, with applications to linear algebra problems[END_REF] to build a new sequence that converges faster towards x. This method works as follows (a more formal definition will be given in Section 3.2): it computes a series of sequences (ε k n ) for k = 1, 2, . . . such that each sequence (ε k n ) for k even converges towards s and the diagonal (d n ) = (ε n 0 ) also converges towards s. This diagonal sequence is the result of the ε-algorithm and is called the accelerated sequence. It converges faster than the original sequence: in only 8 iterates, it reached the fixpoint and stayed constant (see Figure 2, bold lines).

Using the accelerated sequence. When the accelerated sequence reaches the limit (or is sufficiently close to it), we modify the Kleene iteration and directly jump to the limit. Formally, if the limit is (x 1 , x 1 , x 2 , x 2 , x 3 , x 3 ) and if the current Kleene iterate is X p , we construct the abstract element X whose bounds are x 1 , x 1 , . . . and set X p+1 = X p ∪ X and re-start Kleene iteration from X p+1 . In this way, we remain sound (X p ⊆ X p+1 ) and we are very close to the fixpoint, as X ⊆ X p+1 . In this example, Kleene iteration stopped after 2 steps and reached the same fixpoint as the one obtained without widening and acceleration. Figure 3 shows the original Kleene iteration and the modified one, for the infimum of variable x1. Let us recall that the Kleene iteration needed 55 steps to converge, where the modified iteration stops after 18 steps.

Theoretical frameworks

In this section, we briefly recall the basics of abstract interpretation, with an emphasis on the widening operator. Next we present in more details the theory of sequence transformations. Finally, we give our main

Overview of the abstract interpretation theory

Abstract interpretation is a general method to compute over-approximations of program semantics where the two key ideas are:

• Safe abstractions of sets of states thanks to Galois connections. More precisely let C, ⊑ C be the lattice of concrete states and let A, ⊑ A be the lattice of abstract states. A is a safe abstraction of C if there exists a Galois connexion

C, ⊑ C ---→ ← --- γ α A, ⊑ A , i.e.
there exist monotone maps α and γ such

that ∀c ∈ C, ∀a ∈ A, α(c) ⊑ A a ⇔ c ⊑ C γ(a).
• An effective computation method of the abstract semantics with, in general, a widening operator. The semantics of a program is defined as the smallest solution of a recursive system of semantic equations F . Hence, the abstract program semantics is a set of states X of a lattice A, ⊑ A such that X = F (X) where F is monotone. The solution X is iteratively constructed by X i+1 = X i ⊔ F (X i ), starting from X 0 = ⊥. The value ⊥ denotes the smallest element of A and the operation ⊔ denotes the join operation of A. The sequence (X n ) defines an increasing chain of elements of A. This chain may be infinite, so to enforce the convergence of this sequence, we usually substitute the operator ⊔ by a widening operator ∇, see Definition 3.1, that is an over-approximation of ⊔.

Definition 3.1 (Widening operator [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF]) Let A, ⊑ A be a lattice. The map

∇ : A× A → A is a widening operator iff i) ∀v 1 , v 2 ∈ A, v 1 ⊔ v 2 ⊑ A v 1 ∇v 2 . ii) For each increasing chain v 0 ⊑ A • • • ⊑ A v n ⊑ A • • • of A,
the increasing chain defined by s 0 = v 0 and s n = s n-1 ∇v n is stationary:

∃n 0 , ∀n 1 , n 2 , (n 2 > n 1 > n 0 ) ⇒ s n1 = s n2 .
The widening operator plays an important role in static analysis because, thanks to it, we are able to consider infinite state spaces. As a consequence, many abstract domains are associated with a widening operator. For example the classical widening of the interval domain is defined by:

[a, b]∇[c, d] = a if a ≤ c -∞ otherwise , b if b ≥ d +∞ otherwise .
this is to add intermediate steps among a finite set T ; that is the idea behind the widening with thresholds ∇ T . For the interval domain, it is defined [START_REF] Blanchet | A Static Analyzer for Large Safety-Critical Software[END_REF] by:

[a, b]∇T [c, d] = a if a ≤ c max{t ∈ T : t ≤ c} otherwise , b if b ≥ d min{t ∈ T : t ≥ d} otherwise .
While widening with thresholds gives better results, we are facing with the problem to define a priori the set T . Finding relevant values for T is a difficult task for which, to the best of our knowledge, no automatic solution exists.

Acceleration of convergence

We give an overview of the techniques of acceleration of convergence in numerical analysis [START_REF] Brezinski | Extrapolation Methods-Theory and Practice[END_REF]. The goal of convergence acceleration techniques, also named sequence transformations, is to increase the rate of convergence of a sequence. d(xn,ℓ) = 0. This means that (y n ) is asymptotically closer to ℓ than (x n ). An important notion for a sequence transformation T is its kernel K T which is the set of sequences (x n ) for which T (x n ) is ultimately constant. We now present some acceleration methods that we used in our experimentation. For more details, we refer to [START_REF] Brezinski | Extrapolation Methods-Theory and Practice[END_REF].

The Aitken ∆ 2 -method. It is probably the most famous sequence transformation. Given a sequence (x n ) ∈ Ê AE , the accelerated sequence (y n ) is defined by: ∀n ∈ AE, y n = x n -xn+1-xn xn+2-2xn+1+xn . It should be noted that in order to compute y n for some n ∈ AE, three values of (x n ) are required: x n , x n+1 and x n+2 . The kernel K ∆ 2 of this method is the set of all sequences of the form x n = s + a.λ n where s, a and λ are real constants such that a = 0 and λ = 1 (see [START_REF] Brezinski | Generalizations of aitken's process for accelerating the convergence of sequences[END_REF]). The Aitken ∆ 2 -method is an efficient method for accelerating sequences, but it highly suffers from numerical instability when x n , x n+1 and x n+2 are close to each other.

The ε-algorithm. It is often cited as the best general purpose sequence transformation for slowly converging sequences [START_REF] Wynn | The epsilon algorithm and operational formulas of numerical analysis[END_REF]. From a converging sequence (x n ) ∈ Ê AE with limit ℓ, the ε-algorithm builds the following sequences:

(ε -1 n ) : ∀n ∈ AE, ε -1 n = 0, (2)
(ε 0 n ) : ∀n ∈ AE, ε 0 n = xn, (3)
(ε k n ) : ∀k ≥ 1, n ∈ AE, ε k+1 n = ε k-1 n+1 + ε k n+1 -ε k n -1 (4)
The sequence (ε k n ) is called the k-th column, and its construction can be graphically represented as on Figure 4. The even columns (ε 2k n ) (in gray on Figure 4) converge faster to ℓ. The even diagonals (ε 2k n ) k∈AE ) also converge faster to ℓ. In particular, the first diagonal (circled on Figure 4) converges very quickly to ℓ, and it is the accelerated sequence. Let us remark that in order to compute the n-th element of that sequence, 2n elements of (x n ) are required.

Acceleration of vector sequences. Many acceleration methods were designed to handle scalar sequences of real numbers. For almost each of these methods, extensions have been proposed to handle vector sequences (see [START_REF] Graves-Morris | Extrapolation methods for vector sequences[END_REF] for a review of them). The simplest, yet one of the most powerful, of these methods is the vector ε-algorithm (VEA). Given a vector sequence (x n ), the VEA computes a series of vector sequences (ε k n ) using Equations ( 2)-( 4) where the arithmetic operations + andare computed component-wise and the inverse of a vector v is computed as v -1 = v/(v • v), with / being the component-wise division and • the scalar product. The VEA differs from a component-wise application of the (scalar) ε-algorithm as it introduces relations between the components of the vector: the scalar product v • v computes a global information on the vector v which is propagated to all components. Our experiments show that this algorithm works better than a component-wise application of the ε-algorithm. The kernel K ε of the VEA contains all sequences of the form x n+1 = Ax n + B, where A is a constant matrix and B a constant vector [START_REF] Brezinski | A review of vector convergence acceleration methods, with applications to linear algebra problems[END_REF].

Our contribution

In this section, we combine acceleration methods with the abstract fixpoint computation. Our goal is to be as non-intrusive as possible in the classical iterative scheme. In this way, our method can be implemented with minor adaptations in current static analyzers.

Methodology. As seen in Section 3.1, the Kleene iteration for finding the least fixpoint computes with abstract values from some abstract lattice A. In order to use acceleration techniques on the abstract iterates, we need to extract from the abstract elements X n ∈ A a vector of real numbers. Thus, we obtain a sequence of real vectors that we can accelerate, and we quickly reach its limit. We then construct an abstract element X that corresponds to this limit and use it as a candidate for the least fixpoint. This process of transforming an abstract value into a real vector and back is formalized by the notion of extraction and combination functions that are given in Definition 3.2. Definition 3.2 (Extraction and combination.) Let A, ⊑ A be an abstract domain, and let p ∈ AE. The functions Λ A : A → Ê p and Υ A : Ê p → A are called extraction and combination function, respectively, iff for each sequence X n ∈ A AE that order theoretically converges, i.e. ⊔ n∈AE X n = X for some X ∈ A, then the sequence Λ A (X n ) ∈ Ê p AE converges for the usual metric on Ê p , i.e. lim n→∞ Λ A (X n ) = S, and

X ⊑ A Υ A (S).
Intuitively, these functions transpose the convergence of the sequence of iterates into the theory of real sequences, in such a way that the real sequence does not lose any information. Note that the order on Ê p induced by the usual metric is unrelated with the order ⊑ A on A, so the notion of extraction and combination is different from the notion of Galois connection used to compare abstract domains. For the interval domain I = Á v , where v is the number of variables of the program and Á is the set of floating-point intervals, the extraction and the combination functions are defined in Equation. [START_REF] Brezinski | Extrapolation Methods-Theory and Practice[END_REF].

For other domains, these functions must be designed specifically. For example, we believe that such functions can be easily defined for the octagon abstract domain [START_REF] Miné | Weakly Relational Numerical Abstract Domains[END_REF]: the function Λ associates with a difference bound matrix a vector containing all its coefficients. Special care should be taken in the case of infinite coefficients. More generally, we believe that for domains with a pre-defined shape, the functions Λ and Υ can be easily defined. Note that if there is a Galois connection (α I , γ I ) between a domain A

x 0 0 ε 1 0 x 1 ε 2 0 0 ε 1 1 ε 3 0 x 2 ε 2 1 . . . . . . . . . . . .
Arrows depict dependencies: the element at the beginning of the arrow is required to compute the element at the end. For example, ΛI :

ε 2 0 = ε 0 1 + 1 ε 1 1 -ε 1 0 = x1 + 1 ε -1 2 + 1 ε 0 2 -ε 0 1 -ε -1 1 + 1 ε 0 1 -ε 0 0 = x1 + 1 1 x 2 -x 1 -1 x 1 -x 0
I → Ê 2v
(i1, . . . , iv) → i1, i1, . . . , iV , iV ΥI :

Ê 2v → I (x1, x2, . . . , x2v-1, x2v) → [x1, x2], . . . , [x2v-1, x2v] (5)
Accelerated abstract fixpoint computation. We describe the insertion of acceleration methods in the Kleene iteration process in Algorithm 1. We compute in parallel the sequence (X n ) coming from the Kleene's iteration and the accelerated sequence (y n ) computed from an accelerated method. Once the sequence (y n ) seems to converge, that is the distance between two consecutive elements of (y n ) is smaller than a given value δ, we combine the two sequences. That is we compute the upper bound of the two elements of the current iteration. Note that the monotonicity of the computed sequence (X n ) is still guaranteed.

Algorithm 1 Accelerated abstract fixpoint computation 1: repeat 2:

X i := X i-1 ⊔ F (X i-1 )
3:

y i := Accelerate (Λ A (X 0 ), . . . , Λ A (X i )) 4: if ||y i -y i-1 || ≤ δ then 5: X i := X i ⊔ Υ A (y i ) 6: end if 7: until X i ⊑ X i-1
The use of acceleration methods may be seen as an automatic delayed application of the widening with thresholds. Let us remark that we are not guaranteed to terminate in finitely many iterations: we know that asymptotically, the sequence y i from Algorithm 1 gets closer and closer to the fixpoint, but we are not guaranteed that it reaches it. To guarantee termination of the fixpoint computation, we have to use more "radical" widening thresholds, for example after n applications of the accelerated method. So this method cannot be a substitute to widening, but it improves it by reducing the number of parameters (delay and thresholds) that a user must define.

Experimentation

To illustrate our acceleration methods, we used a simple static analyzer1 working on the interval abstract domain that handles C programs without pointers and associated it with our OCaml library of acceleration methods that transform an input sequence (given as a sequence of values) into its accelerated version. The obtained results are presented in the following sections.

Butterworth order 1

To test the acceleration method, we use a first-order Butterworth filter (see Figure 5, left). This filter is designed to have a frequency response which is as flat as mathematically possible in the band-pass and is often used in embedded systems to treat the input signals for a better stability of the program.

The static analysis of this program using the interval abstract domain defines 10 sequences, two for each variable (x1, xn1, y, u, i). These sequences converge toward the smallest fixpoint after a lot of iterations, our acceleration methods allow to obtain the same fixpoint faster. In this example, we accelerate just the upper bound sequences because the lower ones are constant for all the variables. We next present the result obtained with different methods on the variable x1 only, results obtained with other variables are very alike.

The Aitken ∆ 2 -method. In Figure 5, right, with Kleene iteration and without widening, this program converges in 156 iterations, and we get the invariant [0, 20.0084] for x1. With the Aitken ∆ 2 -method, we obtain only in 3 iterations a value very close to 20.0084, but problems of numerical instabilities prevent the stabilization of the program. However the values of the accelerated sequence stay in the interval [20.0082, 20.0086] between the third and the last iteration (see Figure 6, left), which is a good estimate of the convergent point.

The ε-algorithm. In Figure 6, right, we notice a important amelioration in the computation of the fixpoint, thanks to the ε-algorithm. With this method, the fixpoint of the variable x_1 is approximated with a precision of 10 -6 after exactly 8 iterations, while Kleene iteration needed 156 steps. Remark that to obtain 8 elements of the accelerated sequence we need 16 elements from the initial one. We obtain the same results with the vector ε-algorithm.

Butterworth order 2

An order 2 Butterworth filter is given by the following recurrence equation, where x n is a two-dimensional vector, x n = (x On this program, the results obtained using the interval abstract domain are not stable. To address this problem we have used Fluctuat [START_REF] Goubault | Asserting the precision of floating-point computations: a simple abstract interpreter[END_REF], a static analyzer using a specific abstract domain based on affine arithmetic, a more accurate extension of interval arithmetic. It returns the upper and lower bounds of each variables. We applied the vector ε-algorithm on this example with 3 different values of δ (see Algorithm 1): this gives Figure 7. For example, for the variable x 1 and δ = 10 -3 , the over-approximation of the fixpoint is reached after 26 iterations (6 iterations before re-injection and 20 iterations after). Note that we obtain the same fixpoint as with Kleene iteration. We notice that the performance of the Algorithm 1 does not strongly depend of δ. Until now, we use the acceleration just once (unlike in Algorithm 1), a full implementation of it will probably reduce the number of iterations even more.

Related work

Most of the work in abstract interpretation based static analysis concerned the definition of new abstract domains (or improvements of existing ones), and the abstract fixpoint computation remained less studied. Initial work from Cousot and Cousot [START_REF] Cousot | Comparing the Galois connection and widening/narrowing approaches to abstract interpretation[END_REF] discussed various methods to define widening operators. Bourdoncle [START_REF] Bourdoncle | Efficient chaotic iteration strategies with widenings[END_REF] presented different iteration strategies that help reducing the over-approximation introduced by widening. These methods are complementary to our technique: as explained in Section 3.3, acceleration should be done at the same control point as the one chosen for widening, and does not replace standard widening as the termination of the fixpoint computation is not guaranteed. However, acceleration methods greatly improve widening by dynamically and automatically finding good thresholds.

Gopan and Reps in their guided static analysis framework [START_REF] Gopan | Reps. Lookahead Widening[END_REF][START_REF] Gopan | Guided static analysis[END_REF] also used the idea of computing in parallel the main iterates and a guide that shows where the iterates are going. In their work, the precision of the fixpoint computation is increased by computing a pilot value that explores the state space using a restricted version of the iteration function. Once this pilot has stabilized, it is used to accelerate the main iterates; in a sense, this pilot value is very similar to the value y i of Algorithm 1, but we do not modify the iteration function as done in [START_REF] Gopan | Guided static analysis[END_REF].

Maybe the work that is the closest to ours is the use of acceleration techniques in model checking [START_REF] Bardin | FAST: acceleration from theory to practice[END_REF], that have recently been applied to abstract interpretation [START_REF] Gonnord | Combining widening and acceleration in linear relation analysis[END_REF][START_REF] Leroux | Accelerated data-flow analysis[END_REF]. In this framework, the term acceleration is used to describe techniques that try to predict the effect of a loop on an abstract state: the whole loop is then replaced with just one transition that safely and precisely approximates it. These techniques perform very well for sufficiently simple loops working on integer variables, and gives exact results for such cases. Again, this method is complementary to our usage of acceleration: it statically modifies the iteration function by replacing simple loops with just one transition, while our method dynamically predicts the limit of the iterates. We believe that our method is more general, as it can be applied to many kinds of loops and is not restricted to a specific abstract domain (changing the abstract domain only requires changing the Λ A and Υ A functions).

Variable Kleene

Vector ε-algorithm (Before + After) δ = 10 -3 δ = 10 -4 δ = 10 Before: number of iterations to reach the condition on δ. After: the remaining number of Kleene iterations to reach the invariant using the accelerated result.

Conclusion

We presented in this article, a technique to accelerate abstract fixpoint computations using the numerical acceleration methods. This technique consists in building numerical sequences by extracting, at every iteration, supremum and infimum from every variable of the program. We apply to the obtained sequences the various convergence acceleration methods, that allows us to get closer significantly or to reach the fixpoint more quickly than the Kleene iteration. To make sure that the fixpoint returned by the accelerated method is indeed the fixpoint of the abstract semantics, we re-inject it in the static analyzer. This guarantees us the fast stop of the analyzer with a good over-approximation of the fixpoint. The experiments made on a certain number of examples (linear programs) show a good acceleration of the fixpoint computation especially when we use the ε-algorithm, where the number of iterations is divided by four. Let us note that we have assumed in this article that the sequences of iterates and the corresponding vector sequences converge towards a finite limit. In case of diverging sequences, traditional widening can be used as sequence transformation will not perform as well as for converging ones.

For now, we made the experimentation using two separate programs: one that computes the Kleene iterates, and one that accelerates the sequences. The Algorithm 1 is thus still not fully implemented, its automatization is the object of our current work. The use of the interval abstract domain allows to cover just a small set of programs, our future work will also consist in extending this technique to relational domains such as octagons and polyhedra.

Figure 2 :

 2 Figure 2: Sequences extracted from the program of Figure 1 and their accelerated version.

Figure 3 :

 3 Figure3: Infimum value of x 1 . We only display the iterates 5 to 25. At the 15th iteration, the accelerated value is used as a widening with thresholds, and the iteration stops after 18 steps.

  Formally, let D, d be a metric space, i.e. a set D with a distance d : D → Ê + (D will be Ê or Ê p for some p ∈ AE). The set of sequences over D (denoted D AE ) is the set of functions between AE and D. A sequence (x n ) ∈ D AE converges to ℓ iff we have lim n→∞ d(x n , ℓ) = 0. A sequence transformation is a function T : D AE → D AE (T designs a particular acceleration method) such that whenever (x n ) converges to ℓ then (y n ) = T (x n ) also converges to ℓ and lim n→∞ d(yn,ℓ)

Figure 4 :

 4 Figure 4: The ε-table

Figure 5 :Figure 6 :

 56 Figure 5: The Butterworth program (left) and the sequence of supremum of variable x1 (right).

Figure 7 :

 7 Figure 7: Numbers of iterations needed to reach an invariant.

  1 , x 2 ) T :

	xn+1 =	0.9858 -0.009929 0.00929 1	• xn + u •	0.9929 0.004965	, yn+1 =	4.965e -5 0.01	• xn + 2.482e -5 • u

Note that we only consider two consecutive elements to extrapolate the potential fixpoint. The main drawback with this widening is that it may generate too coarse results by going quickly to infinity. A solution of

This analyzer is based on Newspeak, http://penjili.org/newspeak.htmlhttp://penjili.org/newspeak.html, the authors thank especially Sarah Zennou for her technical help.


Jing Dai
email: jing.dai@supelec.fr

Phulpin Yannick
email: yannick.phulpin@supelec.fr

France Supélec

Alain Sarlette
email: alain.sarlette@ulg.ac.be

Damien Ernst
email: dernst@ulg.ac.be

Impact of Delays on a Consensus-based Primary Frequency Control Scheme for AC Systems Connected by a Multi-Terminal HVDC Grid

This paper addresses the problem of sharing primary frequency control reserves among nonsynchronous AC systems connected by a multi-terminal HVDC grid. We focus on a control scheme that modifies the power injections from the different areas into the DC grid based on remote measurements of the other areas' frequencies. This scheme is proposed and applied to a simplified system in a previous work by the authors. The current paper investigates the effects of delays on the control scheme's effectiveness. The study shows that there generally exists a maximum acceptable delay, beyond which the areas' frequency deviations fail to converge to an equilibrium point. This constraint should be taken into account when commissioning such a control scheme.

Introduction

Frequency stability is one of the major concerns for power system operators [START_REF] Kundur | Definition and classification of power system stability[END_REF]. It deals with the power system's ability to return to its nominal frequency after a severe disturbance resulting in a power imbalance. To maintain frequency stability, system operators have developed frequency control schemes, which are usually classified according to the time scale of their actions [START_REF] Rebours | A survey of frequency and voltage control ancillary services -part I: Technical features[END_REF]. The actions corresponding to the shortest time scale are usually referred to as "primary frequency control". It consists of local automatic adjustment, within a few seconds after the disturbance, of the generators' power output based on locally measured frequency variations. The power margins that a generator provides around its scheduled output are named "primary reserves".

In a synchronous AC system, as the average frequency in the time frame of several seconds can be considered identical everywhere, the efforts of the generators participating in primary frequency control sum up within the area. Therefore, larger systems usually experience lower frequency deviations and have lower costs -per MWh -associated to primary frequency control reserves. This has been a significant motivation for interconnecting regional and national systems to create large-scale power systems, such as the UCTE network.

The development of high voltage direct current (HVDC) systems for bulk power transmission over long distances [START_REF] Grünbaum | FACTS and HVDC Light for power system interconnections[END_REF] and underground cable crossings opens new perspectives for interconnecting nonsynchronous areas. In this context, it is generally expected that the power flows through an HVDC system are set at scheduled values, while frequencies of the AC areas remain independent. In a previous paper [START_REF] Dai | Coordinated primary frequency control among nonsynchronous systems connected by a multi-terminal HVDC grid[END_REF], we discuss the possibility of using the fast power-tracking capability of HVDC converters to share primary reserves among non-synchronous areas connected by a DC grid. In the same paper, a control scheme for the HVDC converters is also proposed. It is based on the body of work on consensus problems and relates the problem of sharing primary reserves between the different AC areas to the problem of making the frequency deviations in these areas converge rapidly to the same value.

Both the theoretical study and the simulation results reported in [START_REF] Dai | Coordinated primary frequency control among nonsynchronous systems connected by a multi-terminal HVDC grid[END_REF] show that, under some restrictive assumptions, the scheme can allow a reduction of the requirements in terms of the primary reserves in every AC area. One of these assumptions is that the information on the frequency of one area is instantaneously available to another area of the HVDC system. However, in practice, both measurement and communication introduce delays, which can reach up to a few seconds [START_REF]UCTE operation handbook[END_REF]. As this may challenge the effectiveness of the control strategy, this paper aims to extend the previous findings to more realistic cases by considering the effects of those delays. The analytical part of our study proves a unique equilibrium point for the system following a step change in the load of one AC area. In addition, we provide, under some restrictive assumptions on the power system, conditions under which this equilibrium point is reached. We also report simulation results showing, among others, that delays can indeed cause stability problems.

The paper is organized as follows. Section 2 describes a multi-terminal HVDC system model. Section 3 recalls the control scheme proposed in [START_REF] Dai | Coordinated primary frequency control among nonsynchronous systems connected by a multi-terminal HVDC grid[END_REF] and adds delays to it. Section 4 analyzes stability of a system with such dynamics. Section 5 presents a benchmark system and simulation results.

Multi-terminal HVDC system model

We consider a system with N AC areas, as shown in Fig. 1. The system has three types of components: a DC grid, N non-synchronous AC areas, and N converters that interface the AC areas with the DC grid.

The model aims to reproduce the main characteristics of the frequency of every AC area over a period of several tens of seconds. For notational convenience, the time dependency of the variables is not explicitly mentioned in the equations (i.e., we write x instead of x(t)) where it is not necessary.

AC area

The model of each AC area consists of two components, an aggregated generator and a load, both connected to the AC side of the HVDC converter.

The mechanical dynamics of the generator for area i, ∀i ∈ {1, . . . , N } is described by the equation of motion

2πJ i df i dt = P mi -P ei 2πf i -2πD gi (f i -f nom,i ) (1)
where f i is AC area i's frequency, and f nom,i its nominal value; P mi and P ei are the mechanical power input and the electrical power output of the generator of area i, respectively; and J i and D gi are the moment of inertia and the damping factor of this generator. The power balance within this area requires

P ei = P li + P dc i (2)
where P li is the load demand of area i and P dc i is the power injection from this area into the DC grid.

The AC area's frequency is regulated by the speed governor of the generator, which observes the rotating speed of the shaft and adjusts P mi accordingly, following

T smi dP mi dt = P o mi -P mi - P max mi σ i f i -f nom,i f nom,i (3)
where σ i is the generator droop, T smi the time constant of the servomotor, P max mi the maximum mechanical power available from the turbine, and P o mi the reference value for P mi when f i = f nom,i . In practice, P o mi is refreshed by the secondary frequency controller at a relatively low pace. For instance, its time constant is a few minutes in the UCTE network [START_REF]UCTE operation handbook[END_REF]. Consequently, we assume in this paper that P o mi remains constant. A static load model is used to represent the load

P li = P o li • (1 + D li (f i -f nom,i )) (4)
where P o li is the power demand at f nom,i and D li the frequency sensitivity factor.

DC network

As the electrical time constant of a DC grid is of the order of several milliseconds [START_REF] Kundur | Power System Stability and Control[END_REF], transient dynamics of the DC grid is not considered in our model.

To take into account the general case where there exist nodes that are not connected to any AC area, we suppose that there are in total M ≥ N nodes in the DC grid and that node i is connected to AC area i via converter i, ∀i ∈ {1, . . . , N }. Then, the power transferred from node i to node j within the DC network, denoted by P dc ij , can be expressed as:

P dc ij = V dc i (V dc i -V dc j ) R dc ij (5)
where V dc i and V dc j are the voltages at nodes i and j, respectively, and R dc ij is the resistance between these two nodes. If nodes i and j are not directly connected, R dc ij is considered equal to infinity. Note that there must be either a direct or an indirect connection between any two nodes, otherwise the DC grid would be made of several parts which are not connected to each other.

The power balance at node i satisfies

M j=1 P dc ij = P dc i for i ≤ N , 0 for i > N . (6)
Replacing P dc ij by (5), we can write [START_REF] Kundur | Power System Stability and Control[END_REF] in matrix form as the nonlinear relation

P dc = diag(V dc 1 , . . . , V dc M ) A V dc (7)
where

• P dc is a vector of length M with the first N components equal to P dc 1 , . . . , P dc N and the last M -N components equal to 0. • V dc is a vector containing the DC voltages V dc 1 , . . . , V dc M .

• the components of matrix A are defined as

[A] ij =      - 1 R ij for i = j , j 1 R ij for i = j .

HVDC Converter

Since a converter is capable of tracking a power reference signal with a time constant of several tens of milliseconds [START_REF] Arrillaga | High Voltage Direct Current Transmission[END_REF], its transient dynamics is not considered here.

Conventionally, in an MT-HVDC system, only one of the converters (indexed by k) regulates the DC voltage, while all the others control the real power exchanged between the AC and the DC sides [START_REF] Lescale | Challenges with multi-terminal UHVDC transmissions[END_REF]- [START_REF] Xu | Multi-terminal DC transmission systems for connecting large offshore wind farms[END_REF]. In fact, converter k plays the role of the slack bus that maintains the power balance within the DC grid, and the DC grid can be considered as a real power exchange center between different AC areas. Without loss of generality, we assume an area numbering such that k = N .

Formally, with the notations introduced in Section 2.2, P dc 1 , . . . , P dc N -1 can be used as control variables to influence the frequencies of the different AC areas, whereas P dc N is determined by the DC grid load flow to maintain the power balance within the DC grid.

Coordinated primary frequency controller

This section first recalls the coordinated control scheme proposed in [START_REF] Dai | Coordinated primary frequency control among nonsynchronous systems connected by a multi-terminal HVDC grid[END_REF]. Then, the delays potentially involved in the application of this control scheme are discussed and their influence on the system dynamics is modeled.

Control scheme

The control scheme proposed in [START_REF] Dai | Coordinated primary frequency control among nonsynchronous systems connected by a multi-terminal HVDC grid[END_REF] is distributed in nature. It is composed of N -1 subcontrollers, one for each HVDC converter except converter N which maintains the voltage of the DC grid. The subcontroller assigned to converter i ∈ {1, . . . , N -1} modifies the value of P dc i such that

dP dc i (t) dt =α N j=1 b ij (∆f i (t) -∆f j (t)) + β N j=1 b ij df i (t) dt - df j (t) dt (8)
where

• ∆f i (t) = f i (t) -f nom,i .
• α and β are analogous to integral control gain and proportional control gain, respectively. The influence of their value on the system dynamics will be discussed later in this paper. • b ij is the coefficient representing the communication graph of the system. The value of b ij equals 1 if subcontroller i receives frequency information on area j, and 0 otherwise.

Effects of delays

Equation ( 8) should in principle determine the evolution of P dc i , which represents the power injected by area i into the HVDC grid. However, a subcontroller based on [START_REF] Lescale | Challenges with multi-terminal UHVDC transmissions[END_REF] would lead in a real power system to a variation of P dc i that may differ significantly from the one defined by this equation.

We believe indeed that the sum of the time necessary to measure the frequencies of the AC areas, transmit these values to subcontroller i, compute a reference value for P dc i , and apply it to the converter may be significant. This may in turn lead to an effective variation of P dc i that is delayed with respect to the one predicted by [START_REF] Lescale | Challenges with multi-terminal UHVDC transmissions[END_REF].

To give numerical values for these delays, let us note that it takes at least one period, which is around 20ms (resp. 17ms), to measure a frequency close to 50Hz (resp. 60Hz). Concerning the time necessary to encode, transmit, and decode the frequency information from one AC area to another, it is of the order of several hundreds of milliseconds if not one or two seconds. By way of example, the UCTE does not guarantee delays less than 2 seconds for transmitting information from a substation to a remote SCADA system [START_REF]UCTE operation handbook[END_REF]. As to the time necessary for a converter to effectively inject into the HVDC grid a power setting computed by its subcontroller, it can reach up to tens of milliseconds.

In the following, we will study the properties of the control scheme in the presence of such delays. To simplify the study, we will assume that the overall delays are the same regardless of the subcontroller considered and we denote it by τ .

We will also assume that τ affects the power injected by converter i into the HVDC grid in such a way that the dynamics of P dc i is now given by the following equation:

dP dc i (t) dt =α N j=1 b ij (∆f i (t -τ ) -∆f j (t -τ )) + β N j=1 b ij df i (t -τ ) dt - df j (t -τ ) dt . (9)

System stability

This section reports a theoretical study on the stability properties of the system in the presence of a delay τ . We prove that for a system subjected to a small step change in the load, there exists a unique equilibrium point, at which the frequency deviations are equal in all AC areas. Then, we study the conditions under which the system converges to that equilibrium point, and we provide a Nyquist stability criterion for the special case where all the AC areas have identical parameters.

The theoretical analysis relies on the following assumptions:

Assumption 1: The losses within the DC grid do not vary with time, i.e., N i=1

dP dc i dt = 0 . ( 10
)
Assumption 2: The communication graph that represents the frequency information availability at different subcontrollers has the following properties:

• If the subcontroller of one area has access to the information on another area's frequency, then the subcontroller of this second area also has access to the information on the first area's frequency, i.e., if b ij = 1, then b ji = 1, ∀i = j, i, j ∈ {1, . . . , N -1}. Observe that this property, together with Assumption 1, implies that the time derivative of P dc N also satisfies [START_REF] Jiang | Multiterminal HVDC system in urban areas of large cities[END_REF], where b N i = b iN , ∀i ∈ {1, . . . , N }. • It is constant in time.

Assumption 3: The nonlinear equation ( 1) can be linearized around f i = f nom,i as:

2πJ i df i dt = P mi -P o li -P dc i 2πf nom,i -2πD gi (f i -f nom,i ) .
(11) In the following, the HVDC system is modeled by the linear model, where the dynamics of area i ∈ {1, . . . , N } is defined by (3), [START_REF] Jiang | Multiterminal HVDC system in urban areas of large cities[END_REF], and (11).

Equilibrium point

Proposition 1: Consider that the HVDC system, operating at its nominal equilibrium, is suddenly subjected to a step change in the load of one of its AC areas. Then, under Assumptions 1, 2, and 3, the (linearized) HVDC system has a unique equilibrium point, at which frequency deviations in all AC areas are equal.

Proof: Prior to the step change in the load of one of the AC areas, each AC area is considered in steady state with its frequency regulated at f nom,i . We denote the steady-state values by the variables with a bar overhead. After the step change in the load, the variables start to change. We introduce the following incremental variables:

x i (t) = f i (t) -f nom,i , y i (t) = P mi (t) -Pmi , u i (t) = P dc i (t) -P dc i , v i (t) = P li (t) -P o li .
By introducing these variables, (3), [START_REF] Jiang | Multiterminal HVDC system in urban areas of large cities[END_REF], and (11) become

dx i (t) dt = -a 1i x i (t) + a 2i y i (t) -a 2i u i (t) -a 2i v i (t) , (12)
dy i (t) dt = -a 3i x i (t) -a 4i y i (t) , ( 13
)
du i (t) dt =α N j=1 b ij (x i (t -τ ) -x j (t -τ )) + β N j=1 b ij dx i (t -τ ) dt - dx j (t -τ ) dt (14)
where

a 1i = D gi /J i , a 2i = 1/(4π 2 f nom,i J i ), a 3i = P max mi /(T smi σ i f nom,i
), and a 4i = 1/T smi . Note that a 1i , a 2i , a 3i , and a 4i are all positive constants. Equations ( 12), (13), and ( 14) describe the closedloop dynamics of AC area i, ∀i ∈ {1, . . . , N }, where the state variables are x i (t), y i (t), and u i (t) and the external input is v i (t). Initially, all the state variables are equal to zero, since they are defined as the incremental values with respect to the initial nominal equilibrium. At t 0 , a step change in load occurs in area m such that:

v i (t) = vm for i = m and t > t 0 , 0 otherwise . ( 15
)
We now search for equilibrium points of the system following the step change in the load. Let (x e i , y e i , u e i ) characterize the state of area i at such an equilibrium point. At this point,

dx i (t) dt = dx i (t -τ ) dt = dy i (t) dt = du i (t) dt = 0 . (16)
Thus, (12), (13), and ( 14) become algebraic equations

0 = -a 1i x e i + a 2i y e i -a 2i u e i -a 2i v i (t > t 0 ) , (17) 0 = -a 3i x e i -a 4i y e i , (18)
0 = α N j=1 b ij (x e i -x e j ) . (19)
Equation ( 19) can be written in matrix form for the entire HVDC system. Define the vector x e = [x e 1 , . . . , x e N ] T and let 0 N (resp. 1 N ) denote the column vector of length N with all components equal to 0 (resp. 1). Then, (19) becomes

0 N = αLx e ( 20
)
where L is the Laplacian matrix of the communication graph. It is defined by

[L] ij = -b ij for i = j , j b ij for i = j . ( 21
)
The Laplacian matrix L of a communication graph satisfying Assumption 2 is symmetric positive semidefinite and its only zero eigenvalue corresponds to eigenvector 1 N . Therefore equilibrium requires that the frequency deviations in all AC areas are equal. Let x e be the value of the frequency deviations at the equilibrium point. From ( 17) and (18), we obtain

y e i = - a 3i a 4i x e , (22)
u e i = - a 1i a 2i + a 3i a 4i x e -v i (t > t 0 ) =    -a1i a2i + a3i a4i x e -vm for i = m , -a1i a2i + a3i a4i x e otherwise . ( 23
)
We see in (23) that if x e can be uniquely determined, then the equilibrium point exists and is unique. Assumption 1 implies that N i=1 u i (t) is a constant. The initial conditions yield that N i=1 u i (0) = 0. Thus, N i=1 u e i = 0. From (23), we see that x e is uniquely determined as

x e = -v m • N i=1 a 1i a 4i + a 2i a 3i a 2i a 4i -1 . ( 24
)
Remark 1: The equilibrium point, which is a static quantity, is of course independent of delay τ . However, delays play a crucial role in the system's convergence towards the equilibrium or not, which we study next.

Remark 2: We assumed above that the step change in the load occurs in only one AC area. However, for the general case where v i (t) changes in more than one area and eventually settles at vi different from 0, it is straightforward to extend the above results to reach a similar conclusion on the existence of a unique equilibrium point, with x e given by

x e = - N i=1 vi • N i=1 a 1i a 4i + a 2i a 3i a 2i a 4i -1
. (25)

Stability of the system with identical AC areas

Theoretically proving stability of the control scheme is not an obvious question. We present a result only for the case where all AC areas are assumed identical.

In this subsection we drop AC area index i when referring to the parameters of these areas. We also define the transfer function

h(s) = a 2 (s + a 4 ) (s + a 1 )(s + a 4 ) + a 2 a 3 ( 26
)
where a 1 = D g /J, a 2 = 1/(4π 2 f nom J), a 3 = P max m /(T sm σf nom ), and a 4 = 1/T sm . Proposition 2: Consider that all AC areas of the HVDC system have identical parameters, and that Assumptions 1, 2, and 3 are satisfied. Denote by λ N and λ 2 , respectively the largest and smallest nonzero eigenvalues of the Laplacian associated to the communication graph (see ( 21)). Then the system is stable, and following a step change in the load it asymptotically converges to the unique equilibrium point of Proposition 1, if the net encirclement of any point on the segment (-1/λ 2 , -1/λ N ) by the Nyquist plot of h(s)(α + βs)e -τ s /s is zero.

Proof: Applying the Laplace transform to (12) and (13), we have

x i (s) = a 2 s + a 1 y i (s) - a 2 s + a 1 (u i (s) + v i (s)) , ( 27
)
y i (s) = - a 3 s + a 4 x i (s) . ( 28
)
Eliminating y i (s) from ( 27) yields

x i (s) = -a 2 (s + a 4 ) (s + a 1 )(s + a 4 ) + a 2 a 3 (u i (s) + v i (s)) , (29)
which, written in matrix form, is

x(s) = -h(s)I N (u(s) + v(s)) . ( 30
)
By following the same procedure, the dynamics of P dc i defined by ( 14) can be expressed in the frequency domain as:

u i (s) = ( α s + β)e -τ s N j=1 b ij (x i (s) -x j (s)) , ( 31
)
which can be written in matrix form as:

u(s) = ( α s + β)e -τ s Lx(s) . ( 32
)
By replacing u(s) in ( 30) by (32), we have

x(s) = -s h(s) (sI N + h(s)(α + βs)e -τ s L) -1 v(s) . (33) Define G τ (s) = -s h(s) (sI N + h(s)(α + βs)e -τ s L) -1 . (34)
Then, G τ (s) is the MIMO transfer function between v(s), the load change vector, and x(s), the frequency deviation vector.

The system defined by ( 33) is asymptotically stable if all the poles of its transfer function G τ (s) are on the open left half-plane. Since h(s) is itself a stable transfer function because of the positiveness of a 1 , a 2 , a 3 , and a 4 , we only have to investigate the zeros of Z τ (s) = sI N + h(s)(α + βs)e -τ s L.

Under Assumption 2, the Laplacian L of the communication graph is positive semidefinite and has a single zero eigenvalue. Thus, L = V DV T where V is an orthogonal matrix (containing eigenvectors of L) and D is diagonal (containing eigenvalues 0 =

λ 1 < λ 2 ≤ λ 3 ≤ . . . ≤ λ N ). Now V T Z τ (s)V = sI N +h(s)(α+βs)e -τ s D has the same zeros as Z τ (s).
A single zero at s = 0 is obtained with eigenvector of λ 1 = 0. The latter however cancels with the zero at s = 0 in the numerator1 of G τ (s). To ensure inputoutput stability, Z τ (s) must be positive definite in the subspace spanned by all other eigenvectors (which we denoted by ω k ), for s in the closed right half-plane. This means that

(sI N + h(s)(α + βs)e -τ s L)ω k =sω k + λ k h(s)(α + βs)e -τ s ω k =(s + λ k h(s)(α + βs)e -τ s )ω k =0 N (35)
with k > 1 may not have solutions in the closed right half-plane. The Nyquist criterion says that this holds if the net encirclement of the point (-1/λ k , 0) by the Nyquist plot of h(s)(α+βs)e -τ s /s is zero. Hence the proposition's requirement. Because the whole system's state is observable from the frequency deviation signals, output (i.e., frequency deviation) stability implies stability of the whole state.

The output corresponding to zero initial conditions and a step input

v(t) = 0 for t < 0, v for t > 0,
is then given by

x(s) = G τ (s)v(s) = 1 s G τ (s) v . (36)
As before, we can diagonalize (36) in the basis of eigenvectors of L. From the previous analysis/conditions, components corresponding to λ k , k > 1, have negative poles and therefore, according to linear systems theory, exponentially decay to zero.

The term 1/s in (36) introduced by v(s) is the only term that does not decay away in the output, and in time-domain it corresponds to a step change which represents the shift by x e of all frequency deviations at equilibrium. Remark 3: The above criterion yields that the system defined by (33) is always stable for τ = 0, which is consistent with the theoretical results in [START_REF] Dai | Coordinated primary frequency control among nonsynchronous systems connected by a multi-terminal HVDC grid[END_REF].

To show this, we denote the argument of f (s) by arg(f (s)) and define

J(s) = λ k h(s)(α + βs)/s = λ k a 2 (s + a 4 )(βs + α) s[(s + a 1 )(s + a 4 ) + a 2 a 3 ] . ( 37
)
The argument of J(s) can be calculated as:

arg(J(s)) = arg(s + a 4 ) + arg(βs + α) -90 • -arg((s + a 1 )(s + a 4 ) + a 2 a 3 ) .
With the definition of all the coefficients in J(s), we have arg((s + a 1 )(s + a 4 ) + a 2 a 3 )

< arg((s + a 1 )(s + a 4 )) = arg(s + a 1 ) + arg(s + a 4 ) , from which arg(J(s)) > arg(βs + α) -90 • -arg(s + a 1 )

> -180 • .

On the other hand, it is straightforward to see that arg(J(s)) < 180 • . Thus, the Nyquist plot of J(s) can not intersect with the negative real axis as s grows from j0 to j∞, where j = √ -1. Therefore, the points on the segment (-1/λ 2 , -1/λ N ) are never encircled, which, according to Proposition 2, implies that the system is always stable for τ = 0.

Simulations

To analyze the effects of the delays on the performances of our control scheme, simulations are conducted on an HVDC system with five non-identical non-synchronous areas, which is described in the first part of this section.

Benchmark system

The benchmark system consists of a multi-terminal HVDC grid connecting five non-synchronous areas. The converter of area 5 is chosen to regulate the DC voltage, whose setting is 100kV. The topology of the DC network is represented in Fig. 2. The communication graph coincides with the network topology, i.e., each edge in the figure also represents a bidirectional communication channel between the two areas it connects. The resistance of the DC links are: R 12 = 139Ω, R 15 = 417Ω, R 23 = 278Ω, R 25 = 695Ω, R 34 = 278Ω, and R 45 = 278Ω. In our simulations, we consider that individual AC areas significantly differ from each other, see the parameters in Table 1.

To observe the system's response to a step change in the load, we assume that all the areas operate originally in steady state at their nominal frequency. Then at time t = 2s, a 5% increase of the value of P o l2 (see ( 4)) is modeled.

The continuous-time differential equations ( 3), ( 9), and [START_REF] Ogata | Modern Control Engineering[END_REF] are integrated in this paper using an Euler method with a time-discretization step of 1ms.

Effects of the delays

Simulations reported in [START_REF] Dai | Coordinated primary frequency control among nonsynchronous systems connected by a multi-terminal HVDC grid[END_REF] show that for τ = 0s, the control scheme [START_REF] Lescale | Challenges with multi-terminal UHVDC transmissions[END_REF] drives the frequency deviations of all the areas to the same value. Additionally, when the frequencies are stabilized, the frequency in area 2 is equal to a value which is closer to f nom,2 than when the DC converters are operated with constant power injection.

In contrast, with delays, simulations show that the frequency deviations may fail to converge to each other. In particular, for given values of α and β, there generally exists a maximum acceptable delay beyond which the AC areas' frequencies exhibit oscillations of increasingly large magnitude. For example, when the controller gains are empirically chosen to be 4.44×10 6 , f 1 still converges despite oscillations when τ = 0.35s and fails to converge when τ = 0.37s, as shown in Fig. 3. For comparison, the evolution of f 1 when τ = 0s is also shown in the same figure.

To determine whether the frequency deviations of all the AC areas converge to each other, we define the following criterion, similar to the error band used in the definition of settling time in control theory [START_REF] Ogata | Modern Control Engineering[END_REF].

Let us denote by ∆f e the common value to which the frequency deviations of all AC areas converge when τ = 0s. We classify the system as convergent as long as after t > 22s, i.e., 20 seconds after the step change in the load, all the AC areas' frequency deviations remain within ± 50mHz around ∆f e , i.e.,

|∆f i (t) -∆f e | ≤ 50mHz, ∀i and ∀t > 20s . ( 38
)
We define max as the largest value of the delay for which (38) is satisfied and we search for a relation that may exist between τ max and the controller gains. To ease the analysis, we impose that α = β. We compute τ max for different α = β ∈ [1 × 10 6 , 1 × 10 8 ] by a binary search in τ ∈ [0, 4]s. The points in Fig. 4 represent values of τ max corresponding to different α = β. We can see that by decreasing the values of the controller gains, the harmful oscillations introduced by delays can be curbed. For example, if τ is around two seconds for our system, then we have to decrease the controller gains to a value around 1 × 10 6 to avoid convergence problems. However, as pointed out in [START_REF] Dai | Coordinated primary frequency control among nonsynchronous systems connected by a multi-terminal HVDC grid[END_REF], with lower values of the controller gains, more time is needed for the frequency deviations to converge to similar values.

This phenomenon is illustrated here in the context of a power system with delays by the two sets of curves of Fig. 5 that represent the evolution of the frequencies in the five areas of the system for α = β = 1 × 10 6 and for α = β = 1 × 10 5 . Note that we have also represented in these figures the evolution of f 2 when no control scheme is implemented (i.e., when the power injections into the DC network remain constant).

Conclusions

This paper focuses on a previously proposed control scheme to share primary frequency control reserves among non-synchronous systems connected by a multiterminal HVDC grid. We have studied here the effects of delays on the effectiveness of this control scheme. The study is both analytical and empirical. In particular, we have derived, under some restrictive assumptions on the power system, a stability criterion that may be used to compute the maximum acceptable value for the delay so as to ensure that the control scheme does not lead to stability problems. We have also reported simulation results showing that for delays above a threshold value, the control scheme may cause undamped frequency oscillations. Additionally, as shown by these simulations, these undamped frequency oscillations are more likely to appear when using high values of the controller gains.

As future work, we suggest to extend the theoretical study of the control scheme, notably to the case of non-identical AC areas. We also believe that it would be interesting to test this control scheme on more sophisticated power system benchmarks such as those that would not neglect for example voltage regulation in the AC areas.

Figure 1 .

 1 Figure 1. A multi-terminal HVDC system connecting N AC areas via converters.

•

  The communication graph can not be made of several parts which are not connected to each other, i.e., if b ij = 0, then there must exist some intermediate indices k 1 , . . . , k m such that b ik1 = b k1k2 = . . . = b kmj = 1.

Figure 2 .Figure 3 .

 23 Figure 2. DC grid topology. The circle numbered i represents the point in the DC grid to which converter i is connected. An edge between two circles represents a DC line.

Figure 4 .

 4 Figure 4. Values of τ max for several values of α = β.

5 Figure 5 .

 55 Figure 5. Frequencies of the five AC areas when τ = 2s for α = β = 1 × 10 6 (on the left) and for α = β = 1 × 10 5 (on the right). Both figures also include the evolution of f 2 when the power injections into the DC network remain constant.

Table 1 .

 1 Parameter values for the AC areas.

			Parameter	1	2	Area 3	4	5	Unit
			f nom	50	50	50	50	50	Hz
			P o m P max m J	50 100 2026 6485 6078 2432 2863 80 50 30 80 160 100 60 160	MW MW kg/s 2
			D g	30.5 92.0 88.0 34.5	59.7 kW • s/rad
			σ	0.05 0.10 0.15 0.10 0.075
			T sm	1.5	2.0	2.5	2	1.8	s
			P o l	100	60	40	50	30	MW
			D l	0.01 0.01 0.01 0.01	0.01	s
		10 0
	τ max (s)	10 -1
		10 6 -2 10	10 7 α, β		8 10

This does not correspond to the pole cancellation control. As a matter of fact, the "s" factors in denominator and numerator also cancel for the open-loop system, which is stable. The factors "s" come from rewriting our dynamics, so the pole and zero at s = 0 always cancel exactly.

Acknowledgment

Alain Sarlette is a FRS-FNRS postdoctoral research fellow and Damien Ernst is a FRS-FNRS research fellow. They thank the FRS-FNRS for its financial support. They also thank the financial support of the Belgian Network DYSCO, an Interuniversity Attraction Poles Programme initiated by the Belgian State, Science Policy Office. Alain Sarlette was an invited researcher at the Ecole des Mines de Paris when carrying out this research. The scientific responsibility rests with its authors.