E Wers

H Oudadesse

B Lefeuvre

A Lucas-Girot

J Rocherullé

R Lebullenger

Excess entropy and thermal behavior of Cu-and Ti-doped bioactive glasses

Keywords: Copper, Titanium, Glasses, Entropy, Thermal behavior

Bioactive glasses belong to the ceramic family. They are good materials for implantation due to their excellent capacities to create an intimate bond with bones. Copper is known for its anti-inflammatory, antibacterial, and antifungal properties. Titanium is biocompatible and resistant to corrosion. These chemical elements can be introduced in bioactive glasses to provide a wide variety of uses and to enhance the physiological properties of implanted biomaterials. In this work, bioactive glasses doped with different contents of copper and titanium were synthesized by the melting method. The purpose is to study the effect of doping metal element on the thermal characteristics (T g , T c and T f ). The results revealed that the increase of the content of copper and titanium in the glass matrix decreases the melting temperature and induces an increase of the thermal stability. The excess entropies of pure and doped glasses were calculated. Obtained results highlighted the decrease of the excess entropy with the increase of metal elements contents.

Introduction

Biomaterials are considered like an interesting option to replace bone grafts. Even if the bone is the best biomaterial for filling bone, biomaterials can avoid the risks as the transmission of viruses and are available in unlimited quantity [START_REF] Hench | Bioceramics: from concept to clinic[END_REF][START_REF] Daculsi | Efficacité et performance des substituts osseux pour remplacer les allogreffes et autogreffes[END_REF]. The use of biomaterials was significantly developed for the bone repair. However, the introduction of a new material in the medical field requires a precise study to evaluate the biocompatibility and the physico-chemical properties.

The studied bioactive glasses are used as filling bone biomaterial in orthopedic and maxillofacial surgery. They are composed of several oxides: SiO 2 , Na 2 O, CaO, and P 2 O 5 . Their bioactive character results of their ability to form a hydroxyapatite layer Ca 10 (PO 4 ) 6 (OH) 2 on their surface when they are immersed in a simulated body fluid [START_REF] Wers | Effect of copper and zinc on the bioactivity and cells viability of bioactive glasses[END_REF][START_REF] Oudadesse | Apatite forming ability and cytocompatibility of pure and Zn-doped bioactive glasses[END_REF][START_REF] Dietrich | In vitro bioactivity of melt-derived glass 46S6 doped with magnesium[END_REF]. Moreover, the mineral part of human bones is composed of hydroxyapatite. This formation induces a bone bonding improving the bone growth. Some kind of glasses, in the silicate-phosphate glass compositions, has been studied [START_REF] Kansal | Structural analysis and thermal behavior of diopside-fluorapatite-wollastonite-based glasses and glass-ceramics[END_REF][START_REF] Watts | Influence of magnesia on the structure and properties of bioactive glasses[END_REF][START_REF] Fredholm | Strontium containing bioactive glasses: glass structure and physical properties[END_REF][START_REF] Aina | On the dissolution/reaction of small-grain Bioglass ® 45S5 and F-modified bioactive glasses in artificial saliva (AS)[END_REF][START_REF] Brauer | Fluoride-containing bioactive glasses: fluoride loss during melting and ion release in tris buffer solution[END_REF][START_REF] Goel | Influence of strontium on structure, sintering and biodegradation behavior of CaO-MgO-SrO-SiO2-P2O5-CaF2 glasses[END_REF]. Indeed, Hill et al. [START_REF] Hill | The influence of strontium substitution in fluorapatite glasses and glass-ceramics[END_REF] proved that the adding of strontium in the glass matrix has an impact on the glass transition and crystallization temperatures. Lusvardi et al. [START_REF] Lusvardi | Fluoride-containing bioactive glasses: surface reactivity in simulated body fluids solutions[END_REF][START_REF] Lusvardi | In vitro and in vivo behaviour of zinc-doped phosphosilicate glasses[END_REF] have evaluated the effect of fluoride-and zinc-doped phosphosilicate glasses. Handke, Sitarz, and coauthors [START_REF] Handke | Vibrational spectra of phosphate-silicate biomaterials[END_REF][START_REF] Sitarz | Spectroscopic studies of glassy phospho-silicate materials[END_REF][START_REF] Sitarz | Influence of modifying cations on the structure and texture of silicate-phosphate glasses[END_REF][START_REF] Sitarz | Structure and microstructure of glasses from a NaCaPO 4 -SiO 2 -BPO 4 system[END_REF] have studied the structures of multiples glasses having different chemical compositions. However, there is possible appearance of complications during bone repair such as bacterial infections.

The originality of this study is based on the introduction of Cu and Ti in bioactive glass 46S6 by adding CuO and TiO 2 . Indeed, these chemical elements present interesting physiological properties.

Titanium is mainly used in bio-inert materials such as hip implants made of titanium alloy, screws, and plates to repair bone fractures. It is a chemically biocompatible element which allows having mechanically interesting materials [START_REF] Niinomi | Mechanical properties of biomedical titanium alloys[END_REF].

It is recognized that the use of copper materials in hospitals and in public areas reduces the risk of infections. Copper is a chemical element known for its anti-inflammatory and anti-infectious properties. It is necessary for a lot of biological processes in the human body. Copper is an essential component of the angiogenic response [START_REF] Barbucci | Role of the Hyal-Cu(II) complex on bovine aortic and lymphatic endothelial cells behavior on microstructured surfaces[END_REF]. Angiogenesis is a process in the repair of new biological tissues [START_REF] Hoppe | A review of the biological response to ionic dissolution products from bioactive glasses and glass-ceramics[END_REF]. In vivo copper release has been shown to decrease the risk of ischemia, to induce a vascularised capsule around a cross-linked hyaluronic acid hydrogel [START_REF] Giavaresi | Blood vessel formation after soft-tissue implantation of hyaluronan-based hydrogel supplemented with copper ions[END_REF][START_REF] Erol | Copper-releasing, boron-containing bioactive glass-based scaffolds coated with alginate for bone tissue engineering[END_REF] and to stimulate endothelial cell proliferation [START_REF] Hu | Copper stimulates proliferation of human endothelial cells under culture[END_REF].

The glass is an amorphous system with an unordered structure composed of network forming, modifying, and intermediate oxides. The introduction of metal elements in the glass involves specific changes of thermal characteristics. In a previous work, the effect of zinc has been proven on the melting temperature and the excess entropy. The increase of zinc oxides in the glassy matrix decreases these characteristics [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF]. Copper and titanium oxides are intermediate elements generally used in the industry to give a colored glass. In the glassy matrix, if there are no enough alkaline ions, then the intermediate element will be a network modifier by creating two oxygen bridges. Conversely, if there are enough alkaline ions, then the intermediate element will be a network former [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF][START_REF] Barton | Le verre: science et technologie[END_REF].

Numerous studies have proved the effect of copper and titanium on the thermal behavior and characteristics.

The addition of TiO 2 , into oxides glass-forming system, usually contributes to the stabilization of their structure and to the improvement of their properties: chemical durability, mechanical properties, electrical conductivity [START_REF] Mosner | Thermal properties and stability of lithium titano-phosphate glasses[END_REF].

In the borophosphate glasses, the addition of TiO 2 results in a non-linear increase of glass transition temperature [START_REF] Pospisil | Thermal behaviour and crystallization of titanium-zinc borophosphate glasses[END_REF]. In glasses in the ternary system Na 2 O-TiO 2 -P 2 O 5 , the addition of TiO 2 (from 0 to 5 mol% TiO 2 ) resulted in a non-linear increase of glass transition temperature and dilatation softening temperature. The increase of TiO 2 content involves the increase of the contribution of the surface crystallization mechanism. With increasing TiO 2 content, the temperature of maximum nucleation rate is also gradually shifted from a value close to the glass transition temperature toward the crystallization temperature. The chemical durability of the glasses without titanium oxide is very poor, but with the replacement of Na 2 O or P 2 O 5 by TiO 2 , it increases sharply [START_REF] Mosner | Thermal behaviour and properties of Na 2 O-TiO 2 -P 2 O 5 glasses[END_REF].

In silicate-phosphate glasses, the introduction of copper has an impact on the thermal characteristics.

Indeed, the relationship between the parameters characterizing glass transformation effect and an amount of phosphorous and copper forming the glassy structure showed that an increasing content of phosphorous increases solubility of copper in the structure. Changes of T g were higher in glasses with higher content of P 2 O 5 and CuO [START_REF] Sulowska | Effect of copper addition on glass transition of silicate-phosphate glasses[END_REF]. Moreover, in this kind of glass, the crystallization temperature was found to decrease with the increasing copper content [START_REF] Sulowska | Characterization of thermally induced of crystalline phases in CuO containing silicate-phosphate glasses[END_REF]. Various studies were undertaken to investigate the thermal properties of glasses containing copper. For example, in metallic glasses Cu 50 Zr 43 Al 7 , the glass-forming ability and thermal stability may improve with the adding of Be [START_REF] Lu | Kinetics of non-isothermal crystallization in Cu 50 Zr 43 Al 7 and[END_REF].

Generally, the study of the vitreous state is a approach to the description of the structural disorder. It is the general definition of the entropy. In recent literature on thermal behavior of glasses, several authors believe that the glass formation from a liquid lust be obtained by a loss of entropy at the glass transition temperature [START_REF] Kivelson | Metastable systems in thermodynamics: consequences, role of constraints[END_REF][START_REF] Gupta | On the reality of residual entropies of glasses and disordered crystals[END_REF]. The applicability of thermodynamics to the vitreous state was recognized in the first papers considered now as classic ones [START_REF] Lewis | Zero-point entropy of glasses as physical reality[END_REF]. It was shown that the entropy of a glass is higher than that a crystal of the same composition at T = 0 K [START_REF] Nemilov | Zero-point entropy of glasses as physical reality[END_REF]. Moreover, the temperature during the synthesis of glasses could be an impact on the entropy. Indeed, Tropin et al. [START_REF] Tropin | On the dependence of the properties of glasses on cooling and heating rates: I. Entropy, entropy production, and glass transition temperature[END_REF] have proved that the entropy production in glass transition becomes of increasing importance with an increase of the rates of change of temperature. Moreover, both thermodynamic and statistical-mechanical approaches are shown to confirm the conventional point of view that non-zero residual entropy is frozen-in in the system in the course of the glass transition [START_REF] Tropin | On the dependence of the properties of glasses on cooling and heating rates: I. Entropy, entropy production, and glass transition temperature[END_REF]. The respective value of the residual entropy increases with increasing absolute values of the rate of change of temperature [START_REF] Tropin | On the dependence of the properties of glasses on cooling and heating rates: I. Entropy, entropy production, and glass transition temperature[END_REF]. In a previous work [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF], the entropy calculation was carried out on Zn-doped bioactive glasses. The results have showed that the increase of the amount of Zn in the glass matrix involves a decrease of the excess entropy and an increase of the thermal stability.

Thereby, the purpose of this study is to investigate the effect of the introduction of copper and titanium ions on the thermal characteristics of bioactive glasses used in bone filling. For each chemical composition, the excess entropy was calculated according to the changes of thermal characteristics [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF]. This entropy corresponds to the difference between the melting entropy of crystal and the entropy of glass [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF]. Contrary to crystals, synthesized glasses have entropy different to zero at T = 0 K and versus T f [START_REF] Zarzycki | Les verres et l'état vitreux[END_REF].

Materials and methods

Preparation of bioactive glasses

Eight bioactive glasses had been synthesized from the composition of 46S6 (46 mass% SiO 2 , 24 mass% CaO, 24 mass%, Na 2 O, and 6 mass% P 2 O 5 ). Moreover, this bioactive glass composition 46S6 was studied doped by introduction of different concentrations of copper and titanium: 46S6-xCu and 46S6-xTi (where x = 0, 0.1, 1, 5, and 10 in mass%). TiO 2 and CuO were introduced at the expense of Na 2 O and CaO [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF].

For elaboration of the bioactive glass, sodium metasilicate (Na 2 SiO 3 ), silicon oxide (SiO 2 ), calcium metasilicate (CaSiO 3 ), sodium metaphosphate (Na 3 P 3 O 9 ), copper oxide (CuO), and titanium oxide (TiO 2 ) were weighed and mixed in a polyethylene bottle, for 2 h using a planetary mixer.

The premixed mixtures were melted in platinum crucibles that were placed in an electric furnace. The first rise of temperature rate was 10 °C min -1 and it was hold at 900 °C for 1 h to achieve the decarbonation of all products. The second rise of temperature rate was 20 °C min -1 and it was hold to 1,350 °C for 3 h. The samples were casted in preheated brass molds, in order to form cylinders of 13 mm in diameter, and annealed at 565 °C for 4 h near the glass transition temperature of each glass [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF].

Thermal analysis

Differential thermal analysis (DTA) was used to determine the characteristic temperatures of different bioactive glasses. The DTA principle is based on the detection of whether the phenomenon was exothermic or endothermic phenomenon. The glass transition temperature T g , the crystallization temperature T c and the melting temperature T f have been recorded using a Setaram Labsys 1600TG-DTA/DSC thermal analyzer under N 2 gas atmosphere. Therefore, the onset temperature of crystallization T onset c represented the beginning of the crystallization has been recorded. The bioactive glasses were studied under heating rate of 5 °C min -1 raised from room temperature to 1,400 °C. Moreover, 40 mg of the glass powder was heated in platinum crucible and, at the same time, another empty platinum crucible for use as control. It is necessary to regularly carry out temperature calibration and sensitivity calibration. Thereby, before the thermal analysis, the melting temperature of copper sulfate was checked for calibration. The thermal stability (TS) of bioactive glass has been expressed by the temperature difference between T g and T onset c introduced by Dietzel [START_REF] Dietzel | Glasstruktur und Glaseigenschaften[END_REF]:

TS=T onset c -T g

Results and discussion

Characterizations of bioactive glasses doped with Cu and Ti

The X-ray diffraction patterns were recorded between 5 and 80 (2θ/°) using a Bruker D8 advance diffractometer with Cu Kα radiation. Diffractograms of pure and doped glasses, presented in Fig. 1, revealed a halo of diffraction from 24 to 38 (2θ/°) that was characteristic of an amorphous system. A high content of Cu or Ti introduced in the glass matrix does not affect the amorphous character of these glasses. The infrared spectra were recorded by using Fourier transformer infrared (FTIR) spectrometer Bruker Equinox 55 between 4,000 and 400 cm -1 with a resolution of 2 cm -1 . The infrared spectra of bioactive glasses revealed several characteristic bands (Fig. 2). Fig. 2 FTIR spectra of pure and doped glasses Numerous infrared spectroscopy studies have been carried out about phosphate materials. Jastrzebski et al. [START_REF] Jastrzębski | Infrared spectroscopy of different phosphates structures[END_REF] have shown that non-tetrahedral cations have an influence on the shape of the FTIR spectra and the positions of P-O chemical bonds bands. In silicate glasses, Sitarz [START_REF] Sitarz | The structure of simple silicate glasses in the light of Middle Infrared spectroscopy studies[END_REF] has proven that the presence of alkaline ions in silicate glass structure leads to elimination of Si=O structural defects, which means that the ions are located between the ordered domains. Concerning all chemical compositions, bands characteristic of bend vibrations of Si-O-Si bonds are visible at 479 cm -1 [START_REF] Sitarz | The structure of simple silicate glasses in the light of Middle Infrared spectroscopy studies[END_REF]. Between 1,350 and 1,600 cm -1 , vibrations of C-O bonds are due to the carbon dioxide present during the analysis. At 757 and 924 cm -1 , bands characteristics of tetrahedral vibrations of Si-O-Si bonds are present except for pure glass 46S6 where the band at 757 cm -1 is not present. And a stretch vibration of P=O bonds is visible at 1,050 cm -1 for pure and doped glasses. Thereby, the introduction of doping elements allows creating other chemical bonds.

Changes appear according to the doping element introduced in the glass. At 567 cm -1 , there is a vibration characteristic of bend vibrations of P-O bond for Ti-doped glasses. This bond, in the Ti-doped glasses, could be due to a breaking of P=O bond. Silva et al. [START_REF] Silva | Structural characterization of TiO 2 -P 2 O 5 -CaO glasses by spectroscopy[END_REF] explain that in the glass, with a low content of TiO 2 , the metaphosphate units are dominant and increasing the content; metaphosphate units are replaced by pyrophosphate units. This could also explain the formations P-O bonds for Ti-doped glasses. Indeed, RMN study in calcium phosphate glasses containing titanium has proved that the introduction of TiO 2 in the matrix involves the formation of Ti-O-P bonds [START_REF] Silva | Structural characterization of TiO 2 -P 2 O 5 -CaO glasses by spectroscopy[END_REF]. The increase of TiO 2 could allow the replacement of a P-O-P bond by a Ti-O-P bond and thus increase the intensities of P-O bands at 557 and 1,050 cm -1 .

At 667 and 706 cm -1 , bend vibrations of Si-O-Si bonds are just present on the Cu-doped glasses FTIR spectra. The presence of copper ions in SiO 2 amorphous structure leads to network depolymerization, i.e., breaking of Si-O-Si bridges, leading to an increase in the number of Si-O-non-bridging terminal bonds [START_REF] Sitarz | The structure of simple silicate glasses in the light of Middle Infrared spectroscopy studies[END_REF].

Also, the presence of copper element, whose electronegativity is almost identical to that of silicon, results in additional bands in the FTIR spectra characteristic of Si-O bonds.

Thereby, either copper or titanium, new chemical bonds are formed in the vitreous matrix. Si-O-Si bonds are formed in the Cu-doped glasses while P-O bonds appear in those doped with titanium. The electronegativity of these elements (Table 1) may explain these differences. Copper, with its high electronegativity, will form more new chemical bond than titanium. These data are important to the extent that these formations of chemical bonds may have an influence on the thermal behavior of glasses. The DTA analysis was undertaken in order to measure the characteristic temperatures and calculate the excess entropy.

Table 1 Atomic radius, electronegativity, and fusion point of copper and titanium

Chemical element Atomic radius/pm Electronegativity Fusion point/ °C

Influence of Cu and Ti on the thermal characteristics of bioactive glasses

In bioactive glasses, the ratio CaO/Na 2 O was kept constant and equal to 1. The content of Na 2 O and CaO decreases in favor of the content of CuO or TiO 2 . In a previous work concerning the study of Zn-doped glasses, the variation of the CaO/Na 2 O ratio was involved in the decrease of characteristic temperatures [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF]. Indeed, the introduction of Na 2 O creates two bridge-oxygen in the amorphous system. The two negative charges of oxygen are balanced by the charge of Na + pair forming a neutral electrostatic matrix. Thereby, the network structure is modified and changes the glass thermal properties [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF]. The introduction of CaO does not involve any modifications because the two positive charges of Ca 2+ are balanced and create two tetrahedrons linked by ionic bonds [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF]. Also it is proved that the increase of NaCaPO 4 amount in studied glasses causes a significant increase of their resistance to crystallization [START_REF] Sitarz | Direct crystallization of silicate-phosphate glasses of NaCaPO4-SiO2 system[END_REF].

The impact of Cu and Ti was described.

DTA curves are showed in Fig. 3a for Cu-doped glasses and in Fig. 3b for Ti-doped glasses. Each thermal curve presents 3 peaks characteristics of glass transition, crystallization, and melting temperatures. To determine the start of fusion phenomenon (T onset f) and the start of crystallization phenomenon (T onset c), the tangent method was used, and the thermal characteristics of doped glasses are summarized in Tables 2 and3.

T onset c is an important data to calculate the thermal stability of each glass according to the Dietzel relation. Characteristic temperatures versus the content of Cu or Ti are presented in Figs. 4 and5. It shows that for Cu-doped glasses, a decrease of T g , T c and T f is observed when the content of Cu increases. These results are in accordance with the results obtained by Sulowska et al. [START_REF] Sulowska | Effect of copper addition on glass transition of silicate-phosphate glasses[END_REF]. Indeed, the increasing amount of copper introduced to the amorphous matrix of glasses at the cost of decreasing amount of CaO was the reason of the increasing degree of internal strains relaxation in the glass structure which was accompanied by the decrease of transformation temperature (T g ) [START_REF] Sulowska | Effect of copper addition on glass transition of silicate-phosphate glasses[END_REF]. The same phenomenon is observed for our glasses. Indeed, T g and T c increase (up to 574 and 700 °C respectively), T f decreases (up to 1,085 °C) when the content of Ti increases in the glassy matrix. For T g and T c , these results are found in phosphate glasses [START_REF] Neel | Effect of increasing titanium dioxide content on bulk and surface properties of phosphate-based glasses[END_REF]. Important differences are observed between the low content and the high content of doping element introduced in glasses.

Indeed, between 46S6-0.1Ti and 46S6-10Ti, ΔT g = +26 °C and ΔT c = +89 °C while between 46S6-0.1Cu and 46S6-10Cu, ΔT g = -25 °C and ΔT c = -69 °C. Thereby, the significant result is that Ti involves opposite variations on T g and T c compared to Ag.

And for all glasses, when the amount of metal element increases, the thermal stabilities are higher (up to 129 for 46S6-10Cu and 115 for 46S6-10Ti) as shown in Fig. 6. An important thermal stability reveals a low tendency to crystallization [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF]. Here, Cu-doped glasses present ΔTS = 46 K between 46S6-0.1Cu and 46S6-10Cu. For Ti-doped glasses, ΔTS = 38 K between 46S6-0.1Ti and 46S6-10Ti. Szumera et al. [START_REF] Szumera | Thermal study of Mn-containing silicate-phosphate glasses[END_REF] observed the same influence when MnO2 is introduced in glass matrix of silicate-phosphate glasses. The increase of the amount of MnO2 involves a decrease of the glass transition and an increase of the thermal stability. Finally, the introduction of Cu or Ti in glasses involves a decrease of melting temperature: linear decrease for Cu-doped glasses and non-linear decrease for Ti-doped glasses. For Cu-doped glasses, T c and T g also decrease while for Ti-doped glasses, T g and T c increase. The doping elements introduced in high amount increase the thermal stability of glasses. When the doping metallic oxides are introduced in the amorphous matrix, the covalent Si-O-Si chemical bonds were broken to create metal-oxygen ionic bonds. As covalent bonds are stronger energetically than metal-oxygen bonds formed, changes in the thermal behavior of glasses are observed. Also, as in the works of Sitarz et al. [START_REF] Sitarz | Influence of modifiers and glass-forming ions on the crystallization of glasses of the NaCaPO 4 -SiO 2 system[END_REF] about MgO in phosphosilicate glasses, it would be possible that calcium ions Ca 2+ are replaced by copper ions Cu 2+ causing the reduction of thermal resistance of glasses which would explain the decrease of glass transition, crystallization, and melting temperatures. The changing character of domains structure may be the reason of different thermal behaviors [START_REF] Szumera | Spectroscopic and thermal studies of silicate-phosphate glass[END_REF][START_REF] Szumera | Effect of molybdenum addition on the thermal properties of silicate-phosphate glasses[END_REF]. The chemical characteristics of Ti and Cu (Table 1) may also explain the changes in the thermal behavior. Indeed, the reduction of melting temperatures for Cu-doped glasses is rapid compared to those Ti-doped glasses which start to decrease as from 5 mass% of Ti introduced. The melting points of the introduced metal elements may explain this phenomenon. Copper and titanium have respectively melting points of 1,083 and 1,660 °C that is why the Ti-doped glasses present a slowly decrease of their melting temperature and present a non-linear evolution compared to Cu-doped glasses. Moreover, titanium presents a high electronegativity compared to copper.

Thereby, the Ti-O bonds are stronger than Cu-O bonds. The kinetic of decrease of the melting temperatures could be explained by these electronegativities. A metal having a high electronegativity could involve metaloxygen bonds more important and have an impact on the kinetic of decrease of the melting temperatures.

These data are important in the understanding of thermal behavior of glasses. The introduction of metal elements Cu and Ti have a high impact on the melting temperature and a less impact on the crystallization and glass transition temperatures. In a previous work [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF], the effect of Zn was studied in order to elaborate a protocol of synthesis of a porous biomaterial. This present study can bring information about the behavior of glasses during in vitro assays and their capacity to exchange ionic elements to form hydroxyapatite crystals.

Therefore, we could adapt the use of Cu and Ti in the porous biomaterial to improve its physiological properties and to use it into broader medical applications.

Excess entropy of doped bioactive glasses

The aim of this part is to calculate the excess entropy of each glass and to evaluate the effect of the content of doping metal elements on this value. The entropy at 0 K was calculated. In a previous work [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF], the excess entropy formula was established at a temperature T < T f :

Obtained results show that values obtained for doped bioactive glasses are not equal to zero. However, the excess entropy of pure glass is 79 J K -1 [START_REF] Wers | Thermal behaviour and excess entropy of bioactive glasses and Zn-doped glasses[END_REF]. The obtained excess entropy is from 37 up to 14 J K -1 for Cudoped glasses and from 73 up to 5 J K -1 for Ti-doped glasses. The presence of metal elements Cu or Ti in the amorphous matrix reduces this excess entropy which varies depending on the content of metal element in the glass matrix as shown in Fig. 7. Moreover, up to 3 mass/% of metal element introduced, Cu allows decreasing faster the excess entropy than Ti. Beyond 3 mass% added, there is an inversion and the excess entropy decreases faster with titanium.

Fig. 7 Excess entropy function of the content of doping metal element

The results are in agreement with the Zarzycki theory. We conclude that the third principle of thermodynamic is not applicable to bioactive glasses elaborated in the quaternary system SiO 2 -CaO-Na 2 O-P 2 O 5 . This result will contribute on the comprehension of the changes of the kinetic of bioactivity of doped bioactive glasses compared to pure glass .

The entropy of a crystal is zero at 0 K, i.e., there is more thermal agitation of the atoms within the material and they became motionless. The decrease of the entropy of glasses with the increase of doping elements means that the disorder is reduced in the glass matrix and the atoms present become increasingly immobile. This decrease of the disorder is related to the size of atoms constituting the glass matrix and their electronegativity (property explaining the ability of an atom to attract another and the formation of bonds between atoms).

Titanium is bigger than that of copper (176 pm vs 146 pm). On the contrary, the electronegativity of copper is greater than that of titanium (1.90 vs 1.54). Thereby, the values of excess entropy 46S6-10Cu glasses (14 J K -1 ) and 46S6-10Ti (5 J K -1 ) may be explained by these properties. The size of atoms is a more important characteristic than the electronegativity. Indeed, for Ti-doped glasses, when the titanium content increases, titanium atoms cannot be thermally agitated compared to copper atoms which would reduce the entropy up to a value of 5 J K -1 for 46S6-10Ti. For copper, which has a smaller atomic radius, the entropy decreases by adding of copper, but the entropy of 46S6-10Cu is greater by a factor of 2.8 compared to 46S6-10Ti. The atomic radius being smaller than that of titanium, copper atoms can easily move in the glass matrix hence largest entropy. The capacity of copper created chemical bonds in the glass matrix, due to its important electronegativity, may explain the decrease of entropy in the Cu-doped glasses. The chemical bonds formed tend to order the glass network hence a decrease. These properties of the introduced doping elements are important characteristics to explain the variations of entropy.

Conclusions

The physico-chemical methods of characterization of doped bioactive glasses highlight the amorphous character of all glasses whatever the amount of metal elements introduced in the glass matrix. However, new chemical bonds appear according to the doping element introduced in the glass network.

The introduction of copper or titanium in the glass matrix involves several changes in the thermal characteristics.

Melting temperatures clearly decrease when the amount of these metal elements increases. According to the introduced element, melting temperatures may vary linearly. Changes are observed on the crystallization and glass transition temperatures. The introduction of copper decreases these temperatures, and the introduction of titanium increases them. Moreover, the introduction of metal elements increases the thermal stability of glasses.

The calculated excess entropy shows that the pure glass has more important excess entropy compared to doped glasses. The presence of metal elements in the glass matrix reduces the excess entropy according to the amount of metal introduced. The disorder of the material defined by the entropy decreases with the adding of metal elements in the glass matrix probably due to the formation of metal-oxygen chemical bond that "structure" the glass network.

The present study allows us to understand the thermal behavior of bioactive doped glasses which could have an influence on their chemical reactivity and on their bioactivity. Biomaterials could be adapted in different medical applications by introducing chemically interesting elements having good physiological properties.

Figures Captions

Fig. 1

 1 Fig. 1 Diffractograms of pure and doped glasses

Fig. 3 a

 3 Fig. 3 a Thermal curves of Cu-doped glasses, b thermal curves of Ti-doped glasses

Fig. 4 T

 4 Fig. 4 T c and T g function of the content of doping metal element

Fig. 6

 6 Fig. 6 Thermal stability function of the content of doping metal element

Figure 4 TFigure 5 TFigure 6 Figure 7

 4567 Figure 1 Diffractograms of pure and doped glasses Figure 2 FTIR spectra of pure and doped glasses

Table 2 Thermal

 2

		characteristics of Cu-doped glasses
		46S6-0.1Cu 46S6-1Cu 46S6-5Cu 46S6-10Cu
	L f	57	40	34	22
	T onset c/K 906	893	882	911
	T g /K	823	807	778	752
	TS/K	83	86	104	129
	T f /K	1,486	1,473	1,415	1,358
	ΔS/J K -1 37	28	23	14
	Table 3
	Thermal characteristics of Ti-doped glasses
		46S6-0.1Ti 46S6-1Ti 46S6-5Ti 46S6-10Ti
	L f	110	51	23	9
	T onset c/K 898	915	948	964
	T g /K	821	827	834	849
	TS/K	77	88	114	115
	T f /K	1,475	1,461	1,448	1,358
	ΔS/J K -1 73	36	15	5


Sonia Cafieri

Pierre Hansen

Using mathematical programming to refine heuristic solutions for network clustering

We propose mathematical programming based aproaches to refine graph clustering solutions computed by heuristics. Clustering partitions are refined by applying cluster splitting and a combination of merging and splitting actions. A refinement scheme based on iteratively fixing and releasing integer variables of a mixed-integer quadratic optimization formulation appears to be particularly efficient. Computational experiments show the effectiveness and efficiency of the proposed approaches.

Introduction

Networks, or graphs, provide very useful tools for modelling complex systems [START_REF] Newman | Networks: an introduction[END_REF]. They consist of a set V of vertices associated to the entities under study and a set E of edges each of which joins two vertices and corresponds to relationships among the entities. For instance, in sociology vertices are associated with people and edges with relationships like friendship, communication or collaboration between them. In biology, vertices are associated for instance to proteins and the edges to their interactions. Some topological features of networks are studied to better understand the underlying complex systems, as they may reveal the organizational principles of the system components. The structure of complex systems can in fact be understood by identifying the way the nodes of the corresponding networks are connected to each other. A modular structure characterizes many complex systems, meaning that they contain subgroups of entities sharing some common properties. A topic of particular interest in the study of complex networks is therefore the identification of modules, also called clusters or communities. Given a graph G = (V, E), roughly speaking one seeks subgraphs induced by sets of vertices S i ⊆ V which contain more inner edges (with both vertices in the same subset) than cut edges (with vertices in different subsets). In the last decade the problem of finding clusters in complex networks has been very extensively studied, see Fortunato [START_REF] Fortunato | Community detection in graphs[END_REF] for an in-deep survey.

Many definitions of network modules have been proposed as well as criteria to evaluate partitions of vertices in modules. Maximizing any such criterion over the set of all partitions is a combinatorial optimization problem. The most popular criterion, despite some recent criticism [START_REF] Fortunato | Resolution limit in community detection[END_REF][START_REF] Cafieri | Loops and multiple edges in modularity maximization of networks[END_REF], is the modularity of a subnetwork [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. The modularity of a module is defined as the difference of the fraction of the edges that it contains and the expected number of such edges in a network where edges are distributed at random while keeping the degree distribution of vertices constant, according to the so-called configuration model. Modularity of a partition is the sum of modularities of its clusters. So modularity of a network is a criterion whose maximization provides both the optimal number of clusters and an estimate of the amount of modularity of the network. Numerous heuristics have been proposed for maximizing modularity of a network. They include applications of simulated annealing [START_REF] Guimerà | Functional cartography of complex metabolic networks[END_REF][START_REF] Massen | Identifying communities within energy landscapes[END_REF][START_REF] Medus | Detection of community structures in networks via global optimization[END_REF], mean field annealing [START_REF] Lehmann | Deterministic modularity optimization[END_REF], genetic search [START_REF] Tasgin | Community detection in complex networks using genetic algorithms[END_REF], extremal optimization [START_REF] Duch | Community identification using extremal optimization[END_REF], variable neighborood search [3], spectral clustering [START_REF] Newman | Modularity and community structure in networks[END_REF], linear programming followed by randomized rounding [START_REF] Agarwal | Modularity-maximizing graph communities via mathematical programming[END_REF], dynamical clustering [START_REF] Boccaletti | Detecting complex network modularity by dynamical clustering[END_REF], multilevel partitioning [START_REF] Djidjev | A scalable multilevel algorithm for graph clustering and community structure detection[END_REF], contraction-dilation [START_REF] Mei | Revealing network communities through modularity maximization by a contraction-dilation method[END_REF], divisive [START_REF] Newman | Modularity and community structure in networks[END_REF][START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF] or agglomerative [START_REF] Clauset | Finding community structure in very large networks[END_REF][START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] hierarchical clustering, and several other approaches. Mathematical programming allows us rigorous formulations and solutions for the maximizing modularity optimization problem. Nevertheless, it is rarely used. There are two approaches to use mathematical programming formulations which can be solved to global optimality. Grötschel and Wakabayashi's [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF][START_REF] Grötschel | Facets of the clique partitioning polytope[END_REF] model for clique partitioning can be immediately applied, replacing the original graph by a complete weighted graph. A closed model is used by Brandes et al. [START_REF] Brandes | On modularity clustering[END_REF]. The second approach was proposed by Xu et al. [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF], who express modularity maximization as a mixed-integer quadratic programming problem with a continuous convex relaxation. Column generation can be applied to solve both models [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF]. In these models, modularity is the objective function to be maximized and constraints are used to impose conditions defining a partition of the vertex set. The obtained optimization problems are generally difficult to solve and only small or medium-scale problems can be easily treated . The situation is more favorable when subgraphs of an original graph are handled, as they are more likely to have smaller size (possibly, medium-scale) than the original graph. Given a partition found by a heuristic, one can attempt to refine the result to obtain a new better partition, acting on subnetworks induced by the clusters of the original partition. The purpose of the present paper is to discuss and advance the use of mathematical programming to refine heuristic solutions for network clustering. Two approaches are discussed and compared, one of which is new. The first one was proposed in [START_REF] Cafieri | Improving heuristics for network modularity maximization using an exact algorithm[END_REF] and is based on splitting clusters using an exact algorithm for bipartitioning and merging pairs of clusters. The new one is inspired by the approach in [START_REF] Xu | Module detection in complex networks using integer optimisation[END_REF] and is based on iteratively fixing integer variables and solving the corresponding problem.

The paper is structured as follows. In Sect. 2 we describe the proposed mathematical programming based approaches to refine heuristic partition. In particular, a mixed-integer quadratic model for modularity-maximizing clustering is recalled and the two strategies to refine partitions, that use such a model, are presented. In Sect. 3 a computational analysis and comparison, on a set of instances from the literature, is presented and discussed. Sect. 4 concludes the paper.

Mathematical programming based clustering refinement

Let us consider a partition found by a heuristic for network clustering. It is constituted by subnetworks induced by the clusters found. As a heuristic has been applied, there is no guarantee that the partition given by these subnetworks represents the optimal solution. Thus, one can seek an improved solution by applying a refinement technique. We propose in this section mathematical programming-based refinement techniques, to be employed as post-processing of heuristics for modularity maximization. First, we recall the main elements of a mixed-integer quadratic model for modularity maximization which is used in these refinement techniques.

A MIQP mathematical programming model

Let G = (V, E) be an undirected unweighted graph, with set of vertices V of order n = |V | and set of edges E of size m = |E|. Modularity Q of G can be expressed as the sum of modularities of clusters, each one being a a function of its number of inner edges and of the sum of degrees of its vertices:

Q = ∑ s m s m - D s 2m 2 , (1)
where m s denotes the number of edges in cluster s, and D s denotes the sum of degrees k i of the vertices of cluster s.

In [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] a mixed-integer quadratic formulation is proposed, where (1) is the objective function to be maximized and binary variables are used to identify to which cluster each vertex and each edge belongs. Sets of allocation constraints, and constraints used to express that each vertex belongs to exactly one module, to impose lower and upper bounds on the cardinality of the modules and to break symmetries, fully define the model. In [START_REF] Cafieri | Improving heuristics for network modularity maximization using an exact algorithm[END_REF] this model is specialized to the case of two clusters only, i.e., a bipartition of the graph. Such a model for bipartitioning is recalled below. Notice that it has been also successfully used to build a hierarchical divisive clustering algorithm, see [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF][START_REF] Cafieri | Reformulation of a model for hierarchical divisive graph modularity maximization[END_REF]. First observe that in the case of bipartitioning the objective function ( 1) can be rewritten in a simpler form, expressing the sum of degrees of vertices belonging to one of the two clusters, say D 2 , as a function of the sum of degrees D 1 of vertices belonging to the other one: D 2 = D c -D 1 , where D c denotes the sum of degrees in the cluster c to be bipartitioned. The objective function to split cluster c can then be written as the following quadratic function:

Q c = m 1 + m 2 m - D 1 2 2m 2 - D c 2 4m 2 + D 1 D c 2m 2 . ( 2
)
where m 1 and m 2 are respectively the number of edges inside the two clusters. Decision variables are variables X i, j,s for each edge (v i , v j ) and s = 1, 2, with X i, j,s equal to 1 if the edge (v i , v j ) is inside cluster s and 0 otherwise, and variables Y i,1 for i = 1, 2, . . . n, equal to 1 if the vertex v i is inside cluster 1 and 0 otherwise. Constraints on the problem are allocation constraints, used to impose that any edge (v i , v j ) can belong to cluster s if and only if both of its end vertices i and j also belong to that cluster:

∀(v i , v j ) ∈ E c X i, j,1 ≤ Y i,1 (3) ∀(v i , v j ) ∈ E c X i, j,1 ≤ Y j,1 (4) ∀(v i , v j ) ∈ E c X i, j,2 ≤ 1 -Y i,1 (5) ∀(v i , v j ) ∈ E c X i, j,2 ≤ 1 -Y j,1 (6)
Further constraints express the number of edges of each of the two clusters and the sum of vertex degrees of the first cluster in terms of the decision variables X and Y , and finally integrality constraints are imposed on variables Y . Notice that integrality of variables X is implied by constraints ( 3)-( 6), as well as integrality of D 1 follows by its defining constraint. The following mixed-integer quadratic (MIQP) model, that has a continuous convex relaxation, is finally obtained [START_REF] Cafieri | Improving heuristics for network modularity maximization using an exact algorithm[END_REF]:

(B)                                    max Q c s.t. ∀(v i , v j ) ∈ E c X i, j,1 ≤ Y i,1 ∀(v i , v j ) ∈ E c X i, j,1 ≤ Y j,1 ∀(v i , v j ) ∈ E c X i, j,2 ≤ 1 -Y i,1 ∀(v i , v j ) ∈ E c X i, j,2 ≤ 1 -Y j,1 ∀s ∈ {1, 2} m s = ∑ (v i ,v j )∈E c X i, j,s D 1 = ∑ v i ∈V c k i Y i,1 ∀s ∈ {1, 2} m s ∈ R D 1 ∈ R ∀v i ∈ V c Y i,1 ∈ {0, 1} ∀(v i , v j ) ∈ E c ∀s ∈ {1, 2} X i, j,s ∈ R + 0 .

Splitting and merging clusters

In [START_REF] Cafieri | Improving heuristics for network modularity maximization using an exact algorithm[END_REF] we proposed a refinement technique for clustering results that is built on the mathematical programming formulation (B) recalled above. First, clusters are considered one at a time and the bipartitioning problem (B) is solved exactly, then pairs of clusters are merged and the exact bipartitioning is applied again. More pre-cisely, in a sequence of steps, starting from the original partition obtained applying a heuristic, each cluster is first bipartitioned using an exact algorithm. Notice that (B) is a MIQP with a continuous convex relaxation, that can be solved to global optimality by any standard solver for MIQP problems through the standard branchand-bound method. If the modularity value corresponding to the obtained bipartition is higher than the one of the original cluster, then such original cluster is replaced by the new ones obtained by bipartition, otherwise the original cluster is kept. This sequence of bipartition attempts leads to a new, refined partition. This new partition is furtherly refined by a new sequence of steps, where pairs of clusters, sorted by decreasing number of joining links, are provisionally merged and modularity of the merged cluster is compared to the sum of modularities of the two original clusters. In the case of improvement of the objective function value, the merged cluster is kept at the place of the two original ones. When merging is not beneficial in terms of improvement of the solution, the merged cluster is attempted to be splitted into two parts, according to the procedure applied in the first sequence of refining steps, exactly solving the bipartition problem. The two new clusters are possibly different from the original ones that have been merged, and can potentially correspond to an improved solution.

Fixing integer variables

We now present a novel mathematical programming-based approach to refine heuristic partitions. It is inspired by the methodology proposed by Xu et al. [START_REF] Xu | Module detection in complex networks using integer optimisation[END_REF] for community detection in networks. In [START_REF] Xu | Module detection in complex networks using integer optimisation[END_REF], the authors propose a two-stage procedure, where first a mixed-integer nonlinear problem (similar to that of [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] for a number of clusters generally greater than two, but where the only decision variables are binary variables Y expressing allocation of vertices to modules) is approximately solved to get an initial partition, and then a fixing and releasing scheme is applied.

In this second stage, the authors consider the MIQP model in [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] and solve it, by standard solvers, iteratively fixing a certain number of variables Y to their value 1 and releasing the other variables, that are so free to take a value 1 or 0 depending on the way vertices are re-allocated in the current solution. Fixing integer variables gives a mathematical programming formulation with a reduced number of variables, and so more tractable. We build upon the idea of fixing binary variables, though developing a different approach. Our approach is devised to refine approximate clustering solutions, so we start from the partition provided by a clustering heuristic, that replaces the first stage of the procedure in [START_REF] Xu | Module detection in complex networks using integer optimisation[END_REF]. Then, we attempt to improve the original partition by acting on modules through a new heuristic based on variable fixing. Starting from an assignment of vertices to modules, i.e., from an assignment of 0-1 values to variables Y , we fix n f ix variables to their value 1 and compute a new value for the remaining variables, that is, we re-allocate the corresponding vertices. For each cluster, the vertices that are reallocated are chosen on the basis of their inner degree (the number of neighbors of a vertex inside the cluster), moving first vertices that have a small inner degree and so are likely to have more connections inside a different cluster than the one they are assigned. A given number of (outer) iterations is performed, each one acting on a set Fix, containing variables whose value has to be fixed, and a set Un f ix, containing variables to be released. To avoid using the same sets Fix-Un f ix in successive iterations, random perturbations are applied to these sets.

As acting on the whole graph requires to solve a mixed-integer nonlinear problem that may be quite large even with a number of variables that are fixed, and splitting and merging clusters appears to be an effective strategy for refinements [START_REF] Cafieri | Improving heuristics for network modularity maximization using an exact algorithm[END_REF], we integrate our fixing variables-based strategy in the procedure above based on splitting and merging clusters. To refine a given partition, again we implement the two consecutive steps performing respectively bipartitioning of each cluster and merging mixed to bipartitioning on pairs of clusters. Thus, we consider the MIQP formulation (B), but in place of solving exactly the bipartitioning problem by standard branch-and-bound for MIQP, we apply our fixing variables-based strategy.

Thus, our refinement procedure works as follows. First, each cluster of the original partition is splitted into two sets. To that effect, an initial approximate solution for the bipartition is computed and the above fixing variables-based approach is applied. If the modularity value corresponding to the obtained bipartition is higher than the one of the original cluster, then the original cluster is replaced by the new ones obtained by bipartition, otherwise the original cluster is kept. Once all clusters of the original partition have been examined, the merging-and-splitting procedure is applied. Pairs of clusters, sorted by decreasing number of joining links, are provisionally merged. If merging improves the objective function value, then the merged cluster is kept, otherwise it is splitted into two subsets again applying the fixing variables-based approach.

Computational results

In this section, we apply the proposed clustering refinement techniques to the partitions found by two known and heuristics for graph modularity maximization. The first one was proposed by Noack and Rotta [START_REF] Noack | Multi-level algorithms for modularity clustering[END_REF] and is based on a single-step coarsening with a multi-level refinement. The second one was proposed in 2011 by Cafieri et al. [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF] and is a hierarchical divisive heuristic that is locally optimal in the sense that bipartitions are computed by an exact optimization algorithm.

The first refinement technique (subsection 2.2) is implemented solving the mixedinteger quadratic bipartition problem (B) using CPLEX 12.2 [START_REF]ILOG CPLEX 12.2 User's Manual[END_REF], setting its parameters in such a way that the MIP cutting plane generation is disabled, the branching variable selection strategy is based on reduced pseudo costs, the number of nodes in the Branch-and-Bound tree is limited to 40000, and 1 only thread is used. The fixing variables-based technique (subsection 2.3) is implemented using as a strarting guess an (approximate) affectation of variables provided by CPLEX 12.2 limited to the solution at the root node, and then iterating the fixing variables scheme over 100 iterations. At each iteration, the number of fixed variables is setted to half the cardinality of the current subgraph.

We test the proposed refinement algorithms on datasets in the literature, which correspond to networks modelling various real-life applications. Specifically, we consider a social network of dolphins [START_REF] Lusseau | The bottlenose dolphin community of doubtful sound features a large proportion of long-lasting associations. can geographic isolation explain this unique trait?[END_REF], a network describing interactions among the characters of Hugo's novel Les Misérables [START_REF] Knuth | The Stanford GraphBase: A Platform for Combinatorial Computing[END_REF], a biological network of proteinprotein interactions [START_REF] Dartnell | Selfsimilar community structure in a network of human interactions[END_REF], a network recording co-purchasing of political books on Amazon.com [START_REF] Krebs | [END_REF], a representation of the schedule of games between American college football teams in the Fall of 2000 [START_REF] Girvan | Community structure in social and biological networks[END_REF], a network of connections between US airports [35], a network describing electronic circuits [START_REF] Alon | [END_REF], e-mail interchanges between members of a university [START_REF] Guimerà | Self-similar community structure in a network of human interactions[END_REF], a network giving the topology of the Western States Power Grid of the United States [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF] and authors collaborations [35]. The considered datasets are listed in Table 1 together with their number of vertices n and number of edges m. Solutions have been obtained on a 2.4 GHz Intel Xeon CPU of a computer with 8GB RAM shared by three other similar CPU running Linux.

In Table 2 and Table 3 we report the results of the refinements of clustering results obtained using the Noack and Rotta's [START_REF] Noack | Multi-level algorithms for modularity clustering[END_REF] (NR) heuristic and the Cafieri et al.'s [START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF] (CHL) heuristic respectively. We compare the results of the mathematical programming-based refinements described in Sect. 2, showing the original modularity value computed by the heuristic under consideration (NR or CHL), the intermediate result obtained by cluster splitting only and the final result after sequentially applying the splitting step and the merging step mixed to splitting, for the first refinement technique (split and mrg + spl) (also in [START_REF] Cafieri | Improving heuristics for network modularity maximization using an exact algorithm[END_REF]) and respectively the new one based on fixing variables (split f ix and mrg + spl f ix). We are able to obtain improved results for all the tested cases out of one (political books) refined with the fixing variables technique. Comparing the refined results with optimal modularity maximization solutions, when available in the literature [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF], we remark that in some cases we get the optimal partitions, and in general very good quality solutions. The results obtained applying the two proposed refinements are generally comparable, and often we get the same modularity value (up to 5 decimal digits) in the two cases. When this is not the case, the values coincide up to 2 or 3 decimal digits.

In Table 4 we compare the two proposed approaches in terms of computing time. Very short times are spent in both cases on small-scale networks. For larger networks, it appears that the proposed approach based on fixing integer variables reduces sometimes significantly the time needed to refine the initial partition. This happens, as expected, especially for networks for which exact bipartitioning takes time because of the exploration of a large Branch-and-Bound tree. For example, improving the NR heuristic, time is reduced from 334.72 to 8.96 seconds for the 6-th dataset and from 919.74 to 241.29 seconds for the last dataset, and, improving the CHL heuristic, the reduction is from 454.64 to 16.86 seconds, again for the 6-th dataset. Figure 1 illustrates the clustering of a network for which the optimal modularitymaximizing partition is obtained refining the NR heuristic result.

Conclusions

We proposed mathematical programming based approaches to refine graph clustering solutions. In particular we discussed and compared two approaches, the one in [START_REF] Cafieri | Improving heuristics for network modularity maximization using an exact algorithm[END_REF] based on splitting clusters and a combination of merging and splitting clusters, where bipartitions are computed exactly solving a MIQP problem, and a new one, based on iteratively fixing and releasing integer variables, again integrated in a splitting and merging-splitting scheme. We employ our approach as post-processing of some known heuristics for modularity maximization, obtaining improved solutions and, for some datasets, the optimal partition. The proposed approach based on fix-Table 2 Modularity values corresponding to the partition found by the Noack and Rotta's heuristic [START_REF] Noack | Multi-level algorithms for modularity clustering[END_REF] (Q NR ), by our first approach for refinement after the splitting step only (Q NR split ) and after the merging and splitting step (Q NR mrg+spl ), and by our fixing variables-based approach after the splitting step only (Q NR split f ix ) and after the merging and splitting step (Q NR mrg+spl f ix ). In the last column, the optimal modularity value Q opt is reported, when available in the literature [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF]. (Q CHL ), by our first approach for refinement after the splitting step only (Q CHL split ) and after the merging and splitting step (Q CHL mrg+spl ), and by our fixing variables-based approach after the splitting step only (Q CHL split f ix ) and after the merging and splitting step (Q CHL mrg+spl f ix ). In the last column, the optimal modularity value Q opt is reported, when available in the literature [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF].

dataset Q NR Q NR split Q NR mrg+spl Q NR split f ix Q NR mrg+spl f ix Q opt [
dataset Q CHL Q CHL split Q CHL mrg+spl Q CHL split f ix Q CHL mrg+spl f ix Q opt [

Fig. 1

 1 Fig. 1 Optimal clustering of network football obtained refining the NR heuristic result

Table 1

 1 Datasets in the literature, whith their number of vertices n and number of edges m.

	dataset	n	m
	dolphins	62	159
	les miserables	77	254
	p53 protein	104	226
	political books 105	441
	football	115	613
	usair97	332	2126
	netscience main	379	914
	s838	512	819
	email	1133	5452
	power	4941	6594
	erdos02	6927	11850

Table 3

 3 Modularity values corresponding to the partition found by the Cafieri et al.'s heuristic[START_REF] Cafieri | Locally optimal heuristic for modularity maximization of networks[END_REF]

	2]

Table 4

 4 Computing time (seconds) required by the proposed approaches applied as postprocessing to Noack and Rotta's heuristic (time NR ) and Cafieri et al.'s heuristic (time CHL ). Solutions have been obtained on a 2.4 GHz Intel Xeon CPU of a computer with 8GB RAM shared by three other similar CPU running Linux.

	dataset	time NR mrg+spl	time NR mrg+spl f ix	time CHL mrg+spl	time CHL mrg+spl f ix
	dolphins	0.20	0.39	0.26	0.20
	les miserables	0.67	0.71	0.35	0.30
	p53 protein	1.02	1.23	0.26	0.49
	political books 5.10	1.66	3.41	1.21
	football	3.26	3.16	0.99	0.83
	usair97	334.72	8.96	454.64	16.86
	netscience main	1.38	1.67	0.77	0.85
	s838	1.20	1.40	1.06	1.16
	email	57.80	56.02	-	-
	power	18.62	15.81	17.50	15.42
	erdos02	919.74	241.29	-	-

Acknowledgements The first author has been supported by French National Research Agency (ANR) through grant ANR 12-JS02-009-01 "ATOMIC".


Solenne Moulin

Olivier Pellerin

Loic Toupet
email: loic.toupet@univ-rennes1.fr

Frédéric Paul
email: frederic.paul@univ-rennes1.fr

Formation de palladacycles à six chaînons à partir de précurseurs bisacétate phénantroline Pd(II) et d'isocyanate de phényle Formation of six-membered Palladacycles from Phenanthroline Pd(II) Bisacetate Precursors and Phenylisocyanate

Keywords: Isocyanate, Dimine Pd(II) Complex, Palladacycle, Catlaytic Carbonylation of Nitroaromatics

L'isocyanate de phényle réagit avec des complexes phénanthroline de palladium(II) bisacetate pour former quasi-quantitativement des palladacycles à six chaînons. Dans cette nouvelle réaction, le ligand acétate joue possiblement le rôle d'agent décarbonylant visà-vis de l'isocyanate an formant les palladacycles isolés par un réarrangement intramoléculaire.

Introduction

To date, the reactivity of isocyanates in the coordination sphere of various group-VIII metal centers is still a field only partially explored. [1][2][3] Such reactions can provide a very valuable access to various types of heterocycles and are thus especially attractive with respect to atom-economy considerations. [4][5][6][7] Moreover, they often bring some information about the possible deactivation pathways of catalytic processes involving isocyanate as reactant or product. [8][9][10] In this regard, we have previously shown that diimine-Pd(0) precursors catalytically trimerize arylisocyanates to form triarylisocyanurates (1, Scheme 1), [11] while some inactive palladacycles like 2, 3 or 4a, [9,12] which possibly result from deactivation of the active species, are eventually also formed. Scheme 1. Cyclotrimer and palladacycles formed from diimine-Pd(0) and Pd(II) complexes.

In the process of exploring the scope of this reaction, we have now investigated the reaction of PhNCO with various Pd(II)-phenanthroline precursors, like Pd(N-N)(OAc) 2 complexes (5a-c), where N-N represents 1,10-phenanthroline (o-phen), 3,4,7,8-tetramethylphenanthroline (tmphen) and 2,9-dimethylphenanthroline (dmphen), respectively. In a manner similar to the reaction with Pd(0) precursors, a reaction yielding the six-membered palladacycles 4a-c took place, providing thereby an easy access to these interesting compounds.

Results and discussion

When reacted with one equivalent phenylisocyanate in dichloromethane at 20 °C, the bis-acetato diimine Pd(II) complex 5a[13, 14] forms the known 6-membered palladacycle 4a [9,12] in ca. 30 % yield after two days. In presence of an excess of PhNCO (20 eq.), this reaction becomes nearly quantitative, allowing isolation of ca. 95 % of 4a after the same time. This reaction can also be performed at higher temperature in an aromatic solvent like toluene, 4a being then quantitatively isolated after 16 h (Scheme 2).

Scheme 2. Reaction between 5a and PhNCO.

This reaction can also be performed with the analogue of 5a possessing a tmphen (5b) ligand instead of phen to give the corresponding palladacycle 4b in good yields.

This new palladacycle was characterized by high resolution mass spectrometry (HRMS), infrared and 1 H NMR. In contrast, this reaction is not observed from other Pd(II) precursors such as Pd(o-phen)Cl 2 , Pd(OAc) 2 or from an equimolar Pd(acac) 2 /ophen mixture, pointing at the need to have a carboxylate anion and a diimine ligand simultaneously present in the medium for this reaction to proceed. Notably, we have checked that this reaction takes place with similar yield from a carefully dried complex 5a in anhydrous toluene, suggesting that traces of water are not involved in its mechanism. [START_REF]We have stated that Pd(II) complexes 5a-c often give rise to quite stable solvates with water molecules. A water-free sample of 5a could however be obtained after reacting this complex with acetic anhydride in nitrobenzene (12 h, 80 °C) followed by recrystallization with Et 2 O from nitrobenzene[END_REF] We then turned our attention toward the Pd(II) precursor 5c, presenting the more bulky dmphen ligand in place of o-phen, [14] and with this sterically crowded Pd(II) complex, the reaction takes also place and a nearly complete conversion of 5c into 4c is observed in dichloromethane at ambient temperature after 2 days (Scheme 3). The new palladacycle 4c is the 2,9-dimethylated analogue of 4a. This complex was fully characterized, allowing notably for the observation of the two strong carbonyl  CO stretching modes near 1645 and 1625 cm -1 diagnostic of the three carbonyl groups of this metallacycle. [9] Attempts to perform this reaction at 80 °C in toluene or other aromatic solvents were, however, not successful with partial decomposition of 5c taking place.

Scheme 3. Reaction between 5c and PhNCO.

The steric pressure exerted by the two methyl groups is revealed by the X-ray structure presently obtained for 5cOEt 2 OH 2 (Figure 1). [START_REF]This solvated complex presents nearly identical crystallographic parameters than two different solvates of 5c previously reported[END_REF] It results in a bending of ca. 30° of the phenanthroline core with the mean plane defined by the square planar coordination sphere around the Pd(II) atom. These crystallographic data are in line with those previously obtained for this and the related complex 5a and clearly indicate that 5a-c are mononuclear species. [14,[START_REF] Halligudi | [END_REF] This reaction presents some analogy with the trimerization reaction of PhNCO catalyzed by "(N-N)Pd 0 " species previously studied by some of us. [11] However, the isocyanurate cyclotrimer 1 was never isolated as side-product during any reaction of 5a with PhNCO, suggesting that "(o-phen)Pd 0 " is presently not transiently formed in the medium. Furthermore, the catalytic trimerization of PhNCO was previously shown not to take place with the sterically crowded 2,9-dimethylphenanthroline ligand. When this ligand was used, decomposition of the catalyst into metallic palladium was observed. [11] In contrast, the formation of palladacycle 4c takes readily place at ambient temperatures with the complex 5c. Thus, it is quite likely that the mechanistic pathway for the formation of 4a-c does not involve an elusive "(N-N)Pd 0 " intermediate. ). These steps would be in line with the well established insertion chemistry for this kind of cationic species. [18] It is indeed well known that on such complexes the carboxylate ligands are labile and can be exchanged for other ligands. [9,14] Alternatively, this reaction sequence might be seen as a [Pd(N-N)(OAc)] + -assisted anionic oligomerization of PhNCO by the acetate ion, the latter acting as the nucleophile. This makes sense given the fact that simple acetate salts such as KOAc, in a manner similar to many other inorganic salts, [19] are known to catalyze the cyclotrimerization of arylisocyanates. [20] Moreover, their activity as catalysts in this reaction has been shown to be strongly influenced by the nature of the associated cation. [21] However, this anionicoligomerization process would not ultimately lead to the cyclotrimer 1, but to the palladacycle 4a due to the presence of the second acetate ligand on the nearby Pd(II) center. Indeed, the final steps would involve an intramolecular rearrangement taking place by elimination of acetic anhydride and carbon dioxide to form the metallacyclic core by ring closure. Such an acetate-based intramolecular rearrangement resemble that which was proposed by Moiseev to rationalize the reduction of Pd(II) acetate clusters in presence of CO, [22] except that the "CO" molecule is presently replaced by PhNCO.

Then, depending on the number of inserted isocyanate molecules after ring closure, either 3a-c (n = 1) or 4a-c (n = 2) would be generated in the medium. Since we have previously clearly established that strained palladacycles such as 3a react readily with arylisocyanate to form 4a, [9] this mechanism rationalizes the selective formation of the 6-membered palladacycles 4a-c in presence of excess PhNCO.

[23] Also, based on literature reports, [14,25] the acetate ligand should be more labile in 5c than in 5a. If the first insertion of isocyanate in the Pd-OAc bond in the former complex is rate determining, the reaction might take place as fast with 5c than with 5a, in line with our preliminary observations.

Conclusions

In conclusion we report here a novel and efficient reaction toward the six-membered palladacycles 4a-c from air-stable phenanthroline Pd(II) acetate complexes and phenylisocyanate. While revealing that the particular complexes 5a-c are reactive toward arylisocyanates, this contribution clearly shows that other synthetic routes than these previously evidenced from [(N-N)Pd 0 ] intermediates are possible to obtain palladacycles such as 2-4. Thus, such species might also be directly formed from Pd(II) catalyst precursors when carboxylate ligands are present in the medium. This observation is particularly important in the context of Pd-catalyzed reductive carbonylations of nitroaromatics, because (Pd(II)/diimine/carboxylic acid) catalytic systems are among the most active ones[10] and because the formation of 2-4 and related palladacycles [26] has been proposed to be at the origin of the rapid deactivation of these catalytic systems, [8,9] especially when these reactions are performed in absence of alcohol or amine able to trap the isocyanate formed. [27] The scope and reaction mechanism of this new reaction is currently under investigation.

Experimental

General

The solvents were dried and distilled prior to use and the reaction were performed under inert (Ar) atmospheres in deaerated solvents. Unless otherwise specified, all reagents were purchased from commercial suppliers and used without further purification. Infrared spectra were recorded on a Bruker IFS28 spectrometer (400-4000 cm -1 ).

Reactions of Pd(N-N)(OAc) 2 and PhNCO

Pd(phen){N(C 6 H 5 )C(O)N(C 6 H 5 C(O)N(C 6 H5)} (4a).

To 100 mg of Pd(phen)(OAc) 2 (0.24 mmol) were added 600 mg phenyl isocyanate (20 eq.) dissolved in 15 mL toluene.

The pale-yellow suspension was subsequently heated at 80 °C under stirring. After 16 h, an abundant orange-yellow suspension has formed in the reaction medium. The medium was cooled and 10 mL of ethanol were added to neutralize the excess isocyanate. The suspension is then decanted and washed with several fractions of diethylether yielding 140 mg of yellow-orange solid after subsequent vacuum-drying.

This solid is identified as palladacycle 4a by comparison with authentic samples (95%).

Pd(tmphen){N(C 6 H 5 )C(O)N(C 6 H 5 C(O)N(C 6 H 5 )} (4b).

To 50 mg of Pd(tmphen)(OAc) 2 (0.11 mmol) were added 630 mg phenyl isocyanate (50 eq.) dissolved in 15 mL toluene. The pale-yellow suspension was subsequently heated at 80 °C under stirring. After 16 h, the dark red reaction medium was cooled and 1 mL of EtOH followed by excess diethylether were added. The suspension is then decanted and washed with several fractions of diethylether yielding 39 mg of brown-orange solid after vacuum-drying. This solid is identified as palladacycle 4b (54%

Pd(dmphen){N(C 6 H 5 )C(O)N(C 6 H 5 C(O)N(C 6 H 5 )} (4c).

To 100 mg of Pd(dmphen)(OAc) 2 (0.23 mmol) dissolved in 15 mL of CH 2 Cl 2 were added 550 mg of phenylisocyanate (ca. 20 eq.). The orange solution was stirred 2 days at ambient temperature. The solution is filtrated and diethylether is added to precipitate the title complex as an orange solid. X-Ray Crystallography. Diffraction data frames for 5c were collected on a CCD Saphire 3 Xcalibur apparatus using the graphite monochromatized MoK radiation ( = 0.71073 Å). The cell parameters are obtained with Denzo and Scalepack, [START_REF] Otwinowski | Processing of X-ray diffraction data collected in oscillation mode[END_REF] with 10 frames (psi rotation : 1° per frame). The data collection leads to 4235 independent reflections from which 3338 with I > 2.0(I). [START_REF] Nonius | [END_REF] The structures were solved by direct methods using the SIR97 program, [START_REF] Altomare | [END_REF] and then refined with full-matrix least-square methods based on F 2 (SHELX-97). [START_REF] Sheldrick | SHELX97-2. Program for the refinement of crystal structures[END_REF] The contribution of the disordered solvents (i.

Figure 1 .

 1 Figure 1. ORTEP plot of the palladacycle of 5cOEt 2 OH 2 . Thermal ellipsoids are at the 50% probability level (see Supporting Information for details). Selected bond lengths (Å) and angles (deg): Pd1-O1 1.998(3), Pd1-O2 3.015, Pd1-O3 2.008(3), Pd1-O4 2.963, Pd1-N1 2.021(4), Pd1-N2 2.033(4), C15-O1 1.280(8), C15-O2 1.259(9), C17-O3 1.290(7), C17-O4 1.219(7), N1-Pd1-N2 81.64[START_REF]This solvated complex presents nearly identical crystallographic parameters than two different solvates of 5c previously reported[END_REF], O1-Pd1-O3 85.74[START_REF]This solvated complex presents nearly identical crystallographic parameters than two different solvates of 5c previously reported[END_REF], O1-Pd1-N1 94.58[START_REF]This solvated complex presents nearly identical crystallographic parameters than two different solvates of 5c previously reported[END_REF], O3-Pd1-N2 96.89[START_REF]This solvated complex presents nearly identical crystallographic parameters than two different solvates of 5c previously reported[END_REF].

  Scheme 4. Reaction mechanism proposed for the formation of 4a-c.

CaptionsScheme 1 .Scheme 4 .

 14 Scheme 1. Cyclotrimer and palladacycles formed from diimine-Pd(0) complexes.

Figure 1 .

 1 Figure 1. ORTEP plot of the palladacycle of 5cOEt 2 OH 2 . Thermal ellipsoids are at the 50% probability level (see Supporting Information for details). Selected bond lengths (Å) and angles (deg): Pd1-O1 1.998(3), Pd1-O2 3.015, Pd1-O3 2.008(3), Pd1-O4 2.963, Pd1-N1 2.021(4), Pd1-N2 2.033(4), C15-O1 1.280(8), C15-O2 1.259(9), C17-O3 1.290(7), C17-O4 1.219(7), N1-Pd1-N2 81.64[START_REF]This solvated complex presents nearly identical crystallographic parameters than two different solvates of 5c previously reported[END_REF], O1-Pd1-O3 85.74(16), O1-Pd1-N1 94.58[START_REF]This solvated complex presents nearly identical crystallographic parameters than two different solvates of 5c previously reported[END_REF], O3-Pd1-N2 96.89[START_REF]This solvated complex presents nearly identical crystallographic parameters than two different solvates of 5c previously reported[END_REF].

  Schemes

  NMR (200 MHz, CDCl 3 , , ppm) 8.53 (d, 3 J HH = 8 Hz, 4H, H Ph ), 8.05 (s, 24H, H phen ),

	7.88 (s, 2H, H phen ), 7.53 (t, 3 J HH = 8 Hz, 4H, H Ph ), 7.36 -6.98 (m, 7H, H Ph ), 2.68 (s, 6H,
	CH 3/phen ), 2.20 (s, 6H, CH 3/phen ).

). Color: Brownorange. HRMS (positive ESI, CH 2 Cl 2 /MeOH, m/Z) z 672.1578 [M+H] + , m/z calc'd for [C 36 H 32 N 5 O 2 Pd] 672.1591. FT-IR (KBr, , cm -1 ) 1642 (s, C=O), 1623 (s, C=O). 1 H

  It is collected, washed with several fractions of diethylether and dried in vacuo, yielding 130 mg of the title palladacycle (88 %). Color: HH = 8 Hz, 2H, H phen ), 8.01(d, 3 J HH = 7 Hz, 4H, H Ph ), 7.90 (s, 2H, H phen ), 7.69 (d, 3 J HH = 8 Hz, 2H, H phen ), 7.49 (t, 3 J HH = 8 Hz, 2H, H Ph ), 7.37 (t, 3 J HH = 8 Hz, 4H, H Ph ), 7.26(d, 3 J HH = 8 Hz, 2H, H Ph ), 6.89 (t, 3 J HH = 7 Hz, 4H, H Ph ), 6.69 (t, 3 J HH = 7 Hz, 1 H, H Ph ),

	Orange. Dec. Pt: 180 ± 10 °C. HRMS (positive ESI, CH 2 Cl 2 , m/Z) 644.1291 [M+H] + ,
	m/z calc'd for [C 34 H 28 N 5 O 2 Pd] 644.1278. Anal. Calc'd. for C 34 H 27 N 5 O 2 Pd•H 2 O: C,
	61.68; H, 4.42; N, 10.58; Found: C, 61.67; H, 4.48; N, 10.52. FT-IR (KBr/Nujol, ,
	cm -1 ) 1645 (s, C=O), 1625 (s, C=O). 1 H NMR (200 MHz, CD 2 Cl 2 , , ppm) 8.26 (d,
	2.72 (s, 6H, CH 3/phen ).

3 J

Table of crystallographic parameters for Pd(dmphen)(OAc) 2Table S1 .

 S1 10%). HRMS (positive ESI, CH 2 Cl 2 , m/Z) z 432.0301 [M] + , m/z calc'd for [C 16 H 15 N 2 O 2 Pd] 373.0173. FT-IR (KBr/Nujol, , cm -1 ) 1620 (vs, OAc). 1 H NMR (200 MHz, CD 2 Cl 2 , , ppm) 8.37 (d, 3 J HH = 8 Hz, 2H, H 4 +H 7 ); 7.85 (s, 2H, H 5 +H 6 ), 7.50 (d, 3 J HH = 8 Hz, 2H, H 3 +H 8 ), 2.88 (s, 6H, CH 3/Phen ), 1.93 (s, 6H, OC(O)CH 3 ). Crystal Data, Data Collection, and Refinement Parameters for 5cOEt 2 OH 2 without the solvate(s).

			S5
	Figures for C 18 H 18 N 2 O 4 Pd•H 2 O: C 47.96, H 4.44, N 6.21; Fnd: C 47.51, 4.25, N 6.07. MS number total refl. 4235
	(positive ESI, CH 2 Cl 2 , m/Z) 432 (M + , 3%), 373 (M-OAc, 3%); 286 (M-2OAc, < 1%); number unique refl. 4235
	number obs. refl. [I > 2(I)] R(int) restraints / parameters  = 1/[ 2 (Fo) 2 +(aP) 2 +bP] 209 (dmphen+1, Cpnd formula (where P = [Fo 2 +Fc 2 ]/3) fw final R temp (K) R w	3338 0.000 0/228 a = 0.1 b = 0.0 C 36 H 36 N 4 O 8 Pd 2 0.0583 865.49 120(2) 0.1714	5c
	cryst. syst. R indices (all data)	monoclinic 0.0736
	space group R w (all data) a (Å) Goodness of fit/F 2 (S w )	0.1935 C2/c 18.3528(5) 1.539
	b (Å)  max (e Å -3 ) c (Å)  min (e Å -3 )	18.1223(6) 2.396 Figure 1 13.4071(4) -2.324
	(deg)	90.0
	 (deg)	119.4590(10)
	 (deg)	90.0
	V(Å 3 )	3882.6(2)
	Z	4
	D (calcd) ( g cm -3 )	1.481 (1.796 with solvates)
	crystal size (mm)	0.220.120.08
	F(000)	1744
	abs. coefmm -1 )	0.978
	 range	2.85-27.00
	h k l range	0/23
		0/23
		-17/14

S4

2.

Table of crystallographic parameters for Pd(dmphen)(OAc) 2 p. S4

Acknowledgements

We thank the CNRS for financial support and "Rhone-Poulenc Recherches" for a PhD grant (F.P.) when this work was initiated under the supervision of J.A. Osborn.

Supplementary material

Synthesis and characterization of 5a-d. CCDC-221834 contains the supplementary crystallographic data for this paper. These data can be obtained free of charge from The Cambridge Crystallographic Data Centre via: www.ccdc.cam.ac.uk/data_request.cif.

Supplementary Material

Solenne

Synthesis of Pd(N-N)(OAc) 2 complexes 5a-c

The diimine (N-N) in 10 mL of solvent is added to 250 mg palladium acetate (

Pd(dmphen)(OAc) 2 (5c).

THF was used as solvent with 6.0 equiv. of 2,9dimethylphenanthroline. Yield: 95 %. Dec. Pt: 143 ± 5 °C. Elemental analysis, Calc'd


Sonia Cafieri
email: sonia.cafieri@enac.fr

Maximizing the number of solved aircraft conflicts through velocity regulation *

Keywords: Aircraft conflict avoidance, mixed-integer nonlinear optimization (MINLO), mathematical modelling, Air Traffic Management applications

We propose a model for the maximization of the number of aircraft conflicts that can be solved by performing velocity regulation. The model is mixed-integer as binary variables are used to count solved conflicts and to model alternative choices, while nonlinearities appear in the aircraft separation conditions. The main nonlinearities can however be relaxed by standard reformulations. Numerical results show that the model can be satisfactorly applied at least as a preprocessing step in a conflict avoidance procedure in a given airspace.

Introduction

The growth of air traffic on the world scale leads to an increasing need for automatic decisionsupport tools able to integrate the work of air traffic controllers to guarantee flight safety. In this context, we focus on a crucial problem arising in Air Traffic Management, that of aircraft conflict detection and resolution.

Aircraft are said to be potentially in conflict when their horizontal or altitude distances are less than given standard separation distances (5NM, where 1 NM (nautical mile)= 1852 m, and 1000 ft). Thus, when a loss of separation occurs, aircraft have to be separated. Aircraft conflict avoidance can be performed by different strategies aimed to separate aircraft, including aircraft trajectory (heading angle) changes, flight level changes or aircraft velocity regulation. Corresponding mathematical models can be developed, leading to optimization or optimal control problems. A review is provided in [START_REF] Kuchar | A review of conflict detection and resolution modeling methods[END_REF]. In recent years, mixed-integer linear and nonlinear optimization (MILO, MINLO) have been proposed for aircraft conflict detection and resolution, with interesting results. See for example [START_REF] Alonso-Ayuso | A mixed 0-1 nonlinear optimization model and algorithmic approach for the collision avoidance in ATM: Velocity changes through a time horizon[END_REF], [START_REF] Cafieri | Aircraft deconfliction with speed regulation: new models from mixed-integer optimization[END_REF], [START_REF] Pallottino | Conflict resolution problems for air traffic management systems solved with mixed integer programming[END_REF], [START_REF] Rey | A mixed integer linear model for potential conflict minimization by speed modulations[END_REF]. In previous work [START_REF] Cafieri | Aircraft deconfliction with speed regulation: new models from mixed-integer optimization[END_REF], we proposed MINLO models for conflict avoidance based on velocity regulation, aiming at solving all the conflicts occurring in a given air sector observed during a time horizon (in a tactical flight phase). These models are quite complex and computationally challenging.

In this paper, we propose a mixed-integer nonlinear programming problem for maximizing the number of conflicts that are solved, in a time horizon, when only a velocity regulation is performed. The model can be easily relaxed using standard reformulations. The interest of the proposed model is twofold. First, it allows to easily discriminate between conflicts that can be solved by velocity changes and those that require the application of another separation strategy. Second, it can provide a starting point, and a corresponding feasible solution, for more complex models like those in [START_REF] Cafieri | Aircraft deconfliction with speed regulation: new models from mixed-integer optimization[END_REF], eventually simplifying the branch-and-bound procedure for their solution.

The paper is organized as follows. In Sect. 2 we present the proposed mixed-integer optimization model. In Sect. 3 we discuss the results of some numerical experiments. Sect. 4 concludes the paper.

2.

Model: maximizing the number of solved conflicts

We model conflict avoidance in such a way to achieve aircraft separation by performing a speed change maneuver. This means that a conflict involving a pair of aircraft is solved by aircraft acceleration or deceleration, so that aircraft pass through the points of potential conflict at a different time with respect to what would occur if no maneuvers were performed. There are however a few situations where velocity regulation cannot solve all conflicts of a given aircraft configuration, that corresponds to infeasible optimization problems. In such a case, speed change maneuvers can be performed, leaving potentially some conflicts unsolved and needing the application of another separation maneuver, like heading angle changes.

In the present work we propose an optimization problem where the number of aircraft conflicts that can be solved by speed changes is maximized. The proposed model can then be used as a preprocessing step in a conflict resolution procedure in a target airspace.

Let A be the set of n aircraft. For all i, j ∈ A, let z ij be binary decision variables defined as

z ij = 1 if i and j are separated (no conflict) 0 otherwise
The other decision (continuous) variables of the problem are represented by the aircraft velocities, which are eventually modified with respect to the original ones (that are data of the problem) to solve conflicts:

v min ≤ vi ≤ v max ∀i ∈ A,
where the bounds v min and v max are imposed to allow aircraft only small speed changes, following the idea of subliminal control of velocities suggested in the context of the aeronautic project ERASMUS [START_REF] Bonini | How ERASMUS can support an increase in capacity in 2020[END_REF], such that speeds can vary between -6% and +3% of the original speed. We obtain a mixed-integer model because of the presence of binary as well as continuous variables.

The objective function, to be maximized, is the sum of solved conflicts:

i,j∈A, i<j z ij .
The constraints are given by the integrality constraints on z variables, the bounds on v variables, and the separation constraint on pairs of aircraft. Let us assume that aircraft fly at the same flight level and are identified by 2-dimensional points on a plane. We know their initial position, their trajectory (heading) and their velocity. The aircraft separation between two aircraft i and j at time t is expressed by the condition

||x r ij (t)|| ≥ d, (1)
where d is the minimum required separation distance (usually, 5 NM) and x r ij (t) = x i (t)-x j (t) is a vector representing the relative distance between aircraft i and j. We assume that uniform motion laws apply, so the relative distance of aircraft i and j is expressed as the sum of their relative initial position and the product of their relative speed vr ij by the time:

x r (t) = x r0 ij + vr ij t ∀t,

Table 1 .

 1 Ratio of the aircraft velocities in the optimal solution over the original ones.

	aircraft	vratio
	1	1.00814
	2	1.02809
	3	0.941877
	4	0.981939
	5	0.962551

We proposed a mathematical model for the maximization of the number of aircraft conflicts that can be solved by velocity changes. The model gives a mixed-integer nonlinear optimization problem that can be efficiently solved by standard solvers for MINLO.

* Financial support under grant ANR 12-JS02-009-01 "ATOMIC" is gratefully acknowledged.

that, substituting into (1) and squaring, gives

Notice that the associated equation is an equation of second degree in one unknown t (its graph is a parabola that, as (v r ij ) 2 > 0, has a minimum point and opens upward), that has no solutions if the discriminant

) is negative. The solutions of this equation, if they exist, are the times at which the aircraft are not separated. So, we consider ∆ < 0 as a first condition of separation of aircraft i and j. In the case when this condition is not satisfied, and so aircraft i and j can potentially be in conflict, we look at the form of trajectories. In this work, we assume that trajectories are straight lines intersecting in one point. As per the geometric interpretation of the scalar product, we can look at the sign of the scalar product x r0 ij vr ij to infer if the vectors form an acute or an obtuse angle. In particular, when the scalar product x r0 ij vr ij is negative, then the aircraft are converging, potentially generating a conflict, while they are separated when the product is positive. Finally, we impose aircraft separation imposing that ∆ < 0 or x r0 ij vr ij > 0 for all i, j, i < j. Using again the z binary variables, the two constraints are written as

and, respectively,

Notice that the left hand sides of the two conditions differ only for a square. The same binary variable z ij can be used to model the or condition:

then using an additional variable to account for a separated pair of aircraft when the second condition is satisfied. The nonlinear terms appearing in the constraints come mainly from the products between continuous and binary variables, that can be easily relaxed using the Fortet linearization. This is commonly implemented in the most of the MINLO solvers.

Numerical experiments

We tested our model on instances built placing n aircraft on a circle of a given radius r, in 2-dimensional space, with speed v and a heading angle such that their trajectory is toward the center of the circle (or slightly deviated with respect to such direction). The zone of conflict is around the center of the circle where aircraft are placed, and each aircraft is in conflict with each other. We solve the problem using COUENNE [START_REF] Belotti | Branching and bounds tightening techniques for nonconvex MINLP[END_REF], which implements a spatial Branch-and-Bound based on convex relaxations and provides the global optimal solution.

As an example of solution, let us consider an instance of the conflict avoidance problem with n = 5 aircraft having speed v = 400 NM/h (equal for all aircraft). There are 10 potential conflicts, that are all solved. The ratio of the new speeds over the original ones for the 5 aircraft is shown in Table 1. We see that 2 aircraft are accelerated and 3 of them are decelerated. The speed variation are in the small range [-6%, +3%] around the original velocity for a subliminal control as suggested by ERASMUS. The global optimal solution is obtained in 0.16 seconds.