Cécile Daniel

David Farrusseng

Yves Schuurman
email: yves.schuurman@ircelyon.univ-lyon1.fr

Reaction mechanism of dimethyl carbonate synthesis through isotopic labeling experiments

Keywords: oxygen-18, methanol, carbon dioxide, ceria-zirconia, reaction mechanism

Introduction

Dimethyl carbonate (DMC) is an attractive product because of its low toxicity and fast biodegradability. It is an important synthon for the production of polycarbonates and a methylating agent, but it can also be used as an additive in gasoline, as a solvent in lithium-ion batteries, or as a substitute for VOCs in paints [START_REF] Tan | Review on the synthesis of dimethyl carbonate[END_REF][START_REF] Pacheco | Review of Dimethyl Carbonate (DMC) Manufacture and Its Characteristics as a Fuel Additive[END_REF]. These large-scale applications are bound to increase the global demand for DMC significantly in coming years. DMC was formerly produced on an industrial scale by phosgenation, then more recently by oxidative carbonylation of methanol (Enichem process) or by carbonylation of methyl nitrites (UBE process) [START_REF] Tan | Review on the synthesis of dimethyl carbonate[END_REF][START_REF] Pacheco | Review of Dimethyl Carbonate (DMC) Manufacture and Its Characteristics as a Fuel Additive[END_REF]. These synthesis routes are polluting and toxic. The carboxylation of methanol into DMC appears as a green and sustainable alternative and is produced according to the following equation: However, the thermodynamic equilibrium of this reaction is strongly shifted to the side of the reactants: the conversion of methanol at equilibrium is less than 1% at reaction temperature. But the selective withdrawal of the produced water allows to reach conversions close to 40%. This can be achieved by adding dehydrating reagents to the medium or a process such as catalytic membrane reactors. This requires the understanding of the reaction mechanism and the knowledge of its kinetics. To date, research, mainly devoted to tin-based complexes in homogeneous catalysis and to oxides in heterogeneous catalysis, has not removed the obstacles to promote a high-performance and sustainable catalytic process. The majority of metal oxide studies deal with ceria, zirconia and cerium-zirconium mixed oxides [START_REF] Kohli | Dimethyl Carbonate: Review of Synthesis Routes and Catalysts Used[END_REF][START_REF] Keller | Catalysts, mechanisms and industrial processes for the dimethylcarbonate synthesis[END_REF][START_REF] Tamboli | Catalytic developments in the direct dimethyl carbonate synthesis from carbon dioxide and methanol[END_REF]. A review on DMC synthesis over ceria-based catalysts is given in [START_REF] Tomishige | Catalytic function of CeO2 in nonreductive conversion of CO2 with alcohols[END_REF]. However, the reaction mechanism over ceria-based and related oxide catalysts is still not fully understood.

Bell et al. proposed a complete mechanism over ZrO2, based on infrared and Raman spectroscopic measurements [START_REF] Xie | An in situ Raman study of dimethyl carbonate synthesis from carbon dioxide and methanol over zirconia[END_REF][START_REF] Jung | An in Situ Infrared Study of Dimethyl Carbonate Synthesis from Carbon Dioxide and Methanol over Zirconia[END_REF]. The adsorption of methanol is dissociative and produces water and methoxides species [START_REF] Badri | Use of methanol as an IR molecular probe to study the surface of polycrystalline ceria[END_REF] through the reaction with terminal-or bridged hydroxyl groups. According to [START_REF] Jung | An in Situ Infrared Study of Dimethyl Carbonate Synthesis from Carbon Dioxide and Methanol over Zirconia[END_REF], only the terminal methoxy species are reactive. Next, direct insertion of carbon dioxide into the terminal methoxy species takes place, yielding monomethyl carbonate (MMC). This step does not involve bicarbonate species, which are too stable [START_REF] Jung | An in Situ Infrared Study of Dimethyl Carbonate Synthesis from Carbon Dioxide and Methanol over Zirconia[END_REF][START_REF] Badri | Use of methanol as an IR molecular probe to study the surface of polycrystalline ceria[END_REF]. Less consensus is found for the last step the formation of DMC by the reaction of methanol with MMC. This step is reported to be the rate-limiting step [START_REF] Tomishige | Catalytic properties and structure of zirconia catalysts for direct synthesis of dimethyl carbonate from methanol and carbon dioxide[END_REF][START_REF] Ikeda | Structure of the Active Sites on H3PO4/ZrO2 Catalysts for Dimethyl Carbonate Synthesis from Methanol and Carbon Dioxide[END_REF]. It is commonly proposed that the formation of DMC occurs by the transfer of a methyl species to a MMC species. This implies that the second methanol molecule needs a Lewis acid site to be activated. Ikeda et al. proposed the activation of methanol on a Brønsted acid site on the H3PO4/ZrO2 catalyst [START_REF] Ikeda | Structure of the Active Sites on H3PO4/ZrO2 Catalysts for Dimethyl Carbonate Synthesis from Methanol and Carbon Dioxide[END_REF]. Solid NMR measurements on the catalyst, liquid NMR on the reaction products as well as DRIFT led Aresta et al. [START_REF] Aresta | Influence of Al2O3 on the performance of CeO2 used as catalyst in the direct carboxylation of methanol to dimethylcarbonate and the elucidation of the reaction mechanism[END_REF] to propose the intervention of methanol in the gas phase to form the DMC. DFT studies showed that the step with the lowest energy transition state is the one where two non-adsorbed methanol molecules attack the MMC via an intermolecular step. A study conducted with oxygen-18 labeled methanol (CH3 18 OH) proposed that a methoxide species proceeds via a nucleophilic addition on MMC [START_REF] Honda | Organic carbonate synthesis from CO2 and alcohol over CeO2 with 2-cyanopyridine: Scope and mechanistic studies[END_REF].

Recently, we have reported on zirconia doped ceria synthesized by flame-spray pyrolysis as an active and selective catalyst for the synthesis of DMC [START_REF] Daniel | Surface effect of nano-sized cerium-zirconium oxides for the catalytic conversion of methanol and CO2 into dimethyl carbonate[END_REF]. The highest activities of the CexZr(1-x)O2 catalysts were obtained for molar compositions x ranging between 0.4 and 0.8, but samples prepared by flame-spray pyrolysis were much more active than those prepared by co-precipitation.

The objective of this study is to trace the oxygen atoms during the synthesis reaction of DMC. Do the oxygen atoms in DMC come from the surface and the core of the catalyst (Mars -van Krevelen type mechanism), from carbon dioxide or from methanol? By identifying their origin, steps in the proposed catalytic mechanism can be eliminated or confirmed.

Experimental

Ceria-zirconia as well as pure ceria catalysts were prepared by flame spray pyrolysis as described in [START_REF] Daniel | Surface effect of nano-sized cerium-zirconium oxides for the catalytic conversion of methanol and CO2 into dimethyl carbonate[END_REF]. The flame spray pyrolysis sample, was synthesized following a protocol described in [START_REF] Stark | Flame synthesis of nanocrystalline ceria-zirconia: effect of carrier liquid[END_REF][START_REF] Stark | Flame-made nanocrystalline ceria / zirconia : structural properties and dynamic oxygen exchange capacity[END_REF].

Aqueous metal salt solutions were prepared with desired concentration of metals. The solutions were sprayed into a flame (up to 2200 °C) as a fine mist. High space velocities and heat radiation loss lead to extremely short residence times (milliseconds) with high temperature gradients along the flame axis up to 170 K cm -1 [START_REF] Teoh | Flame spray pyrolysis: An enabling technology for nanoparticles design and fabrication[END_REF]. The metal oxide atoms aggregated into nanoparticles, which were collected on a filter. For this study a sample with the composition Ce0.5Zr0.5O2 was used, referred to as CZ-FSP. Its BET surface area amounted to 99 m 2 /g. Details on its characterization and performance in DMC synthesis are reported in [START_REF] Daniel | Surface effect of nano-sized cerium-zirconium oxides for the catalytic conversion of methanol and CO2 into dimethyl carbonate[END_REF].

To check if the crystal lattice oxygen atoms of CZ-FSP participate during the reaction of DMC synthesis, the oxygen atoms of CZ-FSP where exchanged by oxygen-18 atoms. The isotopic exchange was performed in a fixed bed quartz reactor (4 mm i.d.) with 300 mg of catalyst at 750°C. The isotopic oxygen 18 O2 (Eurisotop) has an isotopic purity of 97.5% vol. The sample was exposed to a flow of 16 O2 and then the flow was switched to a flow of 18 O2 by means of a 4-way valve. The atomic mass units 32, 34 and 36, representative of the oxygen isotopomers, were recorded at the reactor exit as a function of time by a mass spectrometer (Inficon). Once a stable 18 O2 flow was observed, the reactor was cooled down and the sample was recovered. This sample is referred to as 18 CZ-FSP.

The evaluation of catalyst activities was carried out using a method detailed in [START_REF] Daniel | Discovery of very active catalysts for methanol carboxylation into DMC by screening of a large and diverse catalyst library[END_REF]. In the case of the DMC synthesis over the 18 CZ-FSP sample and to minimize the use of CH3 18 OH, a small volume semibatch reactor (4 cm 3 ) was used. The cell loaded with catalyst and CH3OH was placed in an oil bath, regulated at 120°C. The carbon dioxide pressure (C 16 O2 in all cases) was maintained during the 2-hour reaction at 50 bar by a syringe pump (Isco 250D). The measurements were performed with methanol conversions lower than 0.4%, thus far from thermodynamic equilibrium.

At the end of the reaction, the liquid phase was recovered and filtered before being analyzed by a gas chromatograph coupled to a mass spectrometer (GC-MS, Shimadzu QP 2010). Only DMC was detected as product.

The mass spectrum of the unlabeled DMC CH3 16 O (C 16 O) 16 OCH3, with a molar mass (M.M.) of 90 g mol -1 shows the highest fragment 90 amu. DMC CH3 16 O (C 16 O) 18 OCH3 has a molecular fragment at 92 amu instead of 90 amu, and DMC CH3 18 O (C 16 O) 18 OCH3 has a molecular fragment at 94 amu.

Results

Figure 1. 16 O2 to 18 O2 exchange experiment over CZ-FSP at 750°C.

Figure 1 shows the evolution of the oxygen isotopomers as a function of time after the switch from 16 O2 to 18 O2 over CZ-FSP. At first the appearance of 16 O 18 O is observed, followed a little later by 18 O2. This indicates that a hetero-exchange occurs between 18 O2 gas and surface 16 O, in agreement with the results in [START_REF] Madier | Oxygen Mobility in CeO2 and CexZr(1-x)O2 Compounds: Study by CO Transient Oxidation and 18 O/ 16 O Isotopic Exchange[END_REF]. The 16 O 18 O flow goes through a maximum and then levels off to a steady-state production, which corresponds to the 2.5% of 16 O2 present in 18 O2. The 18 O2 flow also reaches steadystate. The exchange was then stopped. Oxygen exchange over 3 monolayers was reported in [START_REF] Madier | Oxygen Mobility in CeO2 and CexZr(1-x)O2 Compounds: Study by CO Transient Oxidation and 18 O/ 16 O Isotopic Exchange[END_REF] for experiments over Ce0.5Zr0.5O2 at 600°C, while other authors report 5 monolayers at 750°C [START_REF] Christou | Suppression of the oxygen storage and release kinetics in Ce0.5Zr0.5O2 induced by P, Ca and Zn chemical poisoning[END_REF].

oxygen-16 atoms due to the hydroxylation of the surface after exposure to air and water contained in methanol. The crystal lattice oxygen atoms are oxygen-18 atoms exchanged over at least 3 monolayers.

The GC-MS analysis of the obtained DMC indicated that it corresponds to CH3 16 O(C 16 O) 16 OCH3. No oxygen-18 atoms were detected in DMC nor in the water produced.

To detect whether 18 O of 18 CZ-FSP participates in the reaction and is inserted into DMC, it is necessary to ensure that the number of oxygen-18 atoms available in 18 CZ-FSP was in large excess with respect to the oxygen involved in the DMC synthesis. We calculated the amount of mobile 18 O in 18 CZ-FSP by assuming that the oxygen-18 atoms in 18 CZ-FSP were available over 3 monolayers with a density of 14 O-atoms nm -2 . This average value is taken from reference [START_REF] Madier | Oxygen Mobility in CeO2 and CexZr(1-x)O2 Compounds: Study by CO Transient Oxidation and 18 O/ 16 O Isotopic Exchange[END_REF] for a ceria-zirconia of equivalent chemical composition and specific surface area. The calculation showed that the amount of 18 O atoms in the 18 CZ-FSP sample during DMC synthesis, 2 mmol, was much larger than the 0.2 mmol of oxygen incorporated in the produced DMC. Thus, there was sufficient oxygen-18 atoms for the production of fully labeled DMC and water.

Another experiment consisted of the reaction between C 16 O2 and CH3 18 This result is consistent with the absence of H2 18 O in the reaction products. Thus, the oxygen atoms of the two OCH3 groups of DMC stem exclusively from CH3 18 OH, not from the surface oxygen-16 hydroxyl groups ( 16 OH).

Discussion

The absence of 18 O in DMC when using an oxygen-18 exchanged ceria-zirconia catalyst implies that the reaction does not proceed according to a Mars van Krevelen type redox mechanism, which consist of a two-step process where oxygen of the catalyst is reacting with one of the reactants or product intermediates, followed by a re-oxidation of the catalyst. The use of CH3 18 OH over the oxygen-16 ceria-zirconia catalyst leads to almost exclusively CH3 18 O (C 16 O) 18 OCH3 (M.M.=94). Honda et al. [START_REF] Honda | Organic carbonate synthesis from CO2 and alcohol over CeO2 with 2-cyanopyridine: Scope and mechanistic studies[END_REF] obtained a ratio of DMC (M.M.=94) to DMC (M.M.=92) of 3.6 for DMC synthesis with CH3 18 OH over pure ceria.

Scheme 1. Sequence of reaction steps for the synthesis of DMC from carbon dioxide and methanol. Adapted from [START_REF] Aresta | Influence of Al2O3 on the performance of CeO2 used as catalyst in the direct carboxylation of methanol to dimethylcarbonate and the elucidation of the reaction mechanism[END_REF].

These results can be interpreted by using the reaction mechanism proposed by Aresta et al. [START_REF] Aresta | Influence of Al2O3 on the performance of CeO2 used as catalyst in the direct carboxylation of methanol to dimethylcarbonate and the elucidation of the reaction mechanism[END_REF],

shown in Scheme 1. In Scheme 1 the 16 O and 18 O atoms have been labeled by a green  and a red *, respectively. Methanol adsorbs dissociatively on a surface hydroxyl group to form a methoxy species and an unlabeled water molecule. Carbon dioxide is inserted into the methoxy species in the next step to form a monomethyl carbonate (MMC) -OCO 18 OCH3 species. Aresta et al. [START_REF] Aresta | Influence of Al2O3 on the performance of CeO2 used as catalyst in the direct carboxylation of methanol to dimethylcarbonate and the elucidation of the reaction mechanism[END_REF] used crosspolarization magic-angle spinning carbon-13 nuclear magnetic resonance spectroscopy to identify the MMC -OCOOCH3 species. They further showed that this species reacts with gas phase methanol, rather than methoxy species to DMC. However, they suggested the reaction path (I) to be the most likely, due to the observed reaction order in methanol that would correspond to the involvement of two methanol-derived species in the rate-controlling step. Moreover, DFT calculations of a similar system showed more favorable energetics for the tri-molecular system than the interaction with a single methanol molecule. On the other hand, Aresta et al. [START_REF] Aresta | Influence of Al2O3 on the performance of CeO2 used as catalyst in the direct carboxylation of methanol to dimethylcarbonate and the elucidation of the reaction mechanism[END_REF] argued that the tri-molecular mechanism may appear entropy forbidden Zirconia doping of ceria is known to increase the oxygen mobility [START_REF] Madier | Oxygen Mobility in CeO2 and CexZr(1-x)O2 Compounds: Study by CO Transient Oxidation and 18 O/ 16 O Isotopic Exchange[END_REF]. It has been reported that the highest oxygen storage capacity of cerium-zirconium mixed oxides, CexZr(1-x)O2, occurs for intermediate compositions 0.4 < x < 0.7 [START_REF] Vidal | Redox behavior of CeO2-ZrO2 mixed oxides: I. Influence of redox treatments on high surface area catalysts[END_REF], just like the maximum activity for DMC synthesis [START_REF] Daniel | Surface effect of nano-sized cerium-zirconium oxides for the catalytic conversion of methanol and CO2 into dimethyl carbonate[END_REF].

However, a direct link between the oxidation state of cerium and the DMC yield has not been clearly established. According to our data, surface or lattice oxygen atoms are not participating to DMC synthesis, ruling out any causal correlation between the oxygen storage capacity and catalyst activity.

A more likely explanation is that zirconia doping of ceria weakens the metal-oxygen bond thus lowering the activation energy of this step.

Conclusions

DMC synthesis from carbon dioxide and methanol have been studied using either an oxygen-18 labeled cerium-zirconium oxide or oxygen-18 labeled methanol. These experiments allowed to shed light on the reaction mechanism. The reaction mechanism consists of three steps: reaction between methanol and surface hydroxyl groups to form water and methoxy species. Reaction between carbon

  OH over oxygen-16 CZ-FSP. The DMC produced was 92% CH3 18 O (C 16 O) 18 OCH3 (M.M.=94). It was identified through the fragment of 94 amu in the GC-MS analysis and by the low intensity of the fragment in 92, representative of CH3 18 O (C 16 O) 16 OCH3. The other larger majority fragments were 31, 33, 47 and 61 amu, also typical of CH3 18 O (C 16 O) 18 OCH3 with a shift of 2 amu's related to the mass of the oxygen-18 isotope.

  [90]. Re-examining the different reaction routes in the light of the results of the labeled methanol experiments, excludes directly the first two routes, (I) & (II) in Scheme 1, for the conversion of the MMC, -OC(O) 18 OCH3 species, as those routes yield H2 18 O and only partially labeled DMC. The route (III), on the other hand, leads to the correctly labeled products DMC and water. DMC is thus formed by a nucleophilic addition of a CH3 18 O methoxide species on MMC -OC(O) 18 OCH3 species. Honda et al. [13] reached a similar conclusion for DMC synthesis over pure ceria. Mixed cerium-zirconium oxides at intermediate compositions are more active than the pure oxides [80]. Zirconia doping of ceria, can promote the catalyst activity either by accelerating the ratecontrolling step or by increasing the number of active sites. Daniel measured apparent activation energies for DMC synthesis over Ce0.5Zr0.5O2 and CeO2, both synthesized by flame spray pyrolysis, of 74 and 86 kJ/mol, respectively [21], shown in Figure 2. This corresponds to approximately a 30-fold increase in activity at 140°C for the zirconia doped sample, which was observed experimentally for samples with similar surface areas [80]. The MMC species have been observed by IR spectroscopy under reaction conditions and seem to be stable intermediates. Step (III) in Scheme 1 is thus likely the rate-controlling step.

Figure 2 .

 2 Figure 2. Arrhenius plot for the synthesis of DMC from methanol and carbon dioxide over Ce0.5Zr0.5O2 and CeO2.

CZ-FSP was then used for the synthesis of DMC with CH3[START_REF] Stark | Flame-made nanocrystalline ceria / zirconia : structural properties and dynamic oxygen exchange capacity[END_REF] OH and C 16 O2. During the reaction, the surface of[START_REF] Daniel | Discovery of very active catalysts for methanol carboxylation into DMC by screening of a large and diverse catalyst library[END_REF] CZ-FSP was composed of the following: the oxygen atoms of the hydroxyl groups were

Acknowledgements

The work leading to these results received funding from the European Union Seventh Framework Program FP7-NMP-2010-Large-4, under grant agreement no. 263007 (acronym CARENA). We thank Johnson-Matthey for providing samples.

dioxide and surface methoxy species to give surface MMC species. Finally, the rate controlling step is a nucleophilic attach of methanol on the MMC species to yield DMC. The oxygen atoms in DMC stem either from methanol or carbon dioxide, but oxygen atoms from cerium-zirconium oxide do not participate in the reaction. The role of zirconia is to weaken the metal-oxygen bond of the MMC species and thus accelerate the reaction by lowering the activation energy of the last step. However, further studies are necessary to confirm this, for example through ab initio calculations.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Jesús Romero González

Fatima Rahmouni

Benjamin Pointard

Andrea Bertoldi

Jérôme Lodewyck

Rodolphe Le Targat

RAZPOUTYNE project towards compact cooling techniques

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The project aims at loading >10 4 atoms into the optical trap in <150 ms.

When interrogating the clock transition the result is judged by two criteria :

• The uncertainty, or how close we believe the corrected measured frequency is to the 'ideal' unperturbed value, and it will be reflected by how well are known the corrections made due to systematical effects shifting this frequency.

• The stability, when the values of repeated measurements are clustered close together and have minimal scatter over time.

The major limitation noise on the stability

Dick effect

Loading as fast as the coherence of the laser will reduce it

Magnetic field along Cooling cross XZ Plane (Y= 0 mm) . B z Field [mT]

Arrangement of permanent magnets ❑ Region 1 and 2 can be use for deceleration (also region 3).

❑ ZS beam contains of 50% σ -and 50% σ + , since the polarization is perpendicular to B field.

❑ Two different circular polarization play a roll:

• σ + light slow down atoms where ∇B>0, and they are capture in the 2D-MOT region.

• σ -light slow down atoms where ∇B<0, and they are capture in the 2D-MOT region after the atomic beam "reflection".

❑ The laser detuning (Δ laser ) able us to select a velocity window, being an offset.

Fundings

ANR-20 CE47-0006 ❑ Along the Y axis:

 The Yb atoms are confinement in a cigarette shape along the center of the MOT.

 A pushing beam will translate the Yb cigarette shape to the 3D-MOT.

❑ Key aspect:

 The same array of magnets create the necessary magnetic field for the Zeeman Slower and the magnetic gradient for the 2D-MOT.

• Zeeman Slower with permanent magnets :

[ 2 ]❑

 2 Lamporesi, Giacomo & Donadello, Simone & Serafini, Simone & Ferrari, G. (2013). Compact high-flux source of cold sodium atoms. The Review of scientific instruments. 84. 063102. 10.1063/1.4808375. . Doppler cooling & Zeeman shift 1D Cool down longitudinal velocity Doppler cooling & Zeeman shift 2D Cool down the radial plane of velocities • 2D-MOT : Distance along optical axis X [mm] = 80 mW • w zs = 12.5 mm • Δ laser = -400 MHz • v cap = 300 m/s • P 2D-MOT = 60 mW • w 2D-MOT = 12.5 mm • Δ 2D-MOT = 0X axis, Zeeman slower beam is counterpropagating to the Yb atoms, who are launched from from the oven. ❑ ZX plane, 4 counter-propagating beam create the 2D-MOT.


domain: info math

text:
Albert Benveniste
email: albert.benveniste@inria.fr

Benoît Caillaud
email: benoit.caillaud@inria.fr

Mathias Malandain
email: mathias.malandain@inria.fr

Joan Thibault
email: joan.thibault@inria.fr

Towards the separate compilation of Modelica: modularity and interfaces for the index reduction of incomplete DAE systems

Keywords: DAE, Modelica, object-oriented modeling, index reduction, structural analysis, linear programming, interface theory, difference bound matrices

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

System modeling tools are key to the engineering of safe and efficient Cyber-Physical Systems (CPS). Although ODE-based languages and tools, such as Simulink (Math-Works, Inc. 1994[START_REF]Modelica, A Unified Object-Oriented Language for Systems Modeling[END_REF], are widely used in industry, there are two main reasons why DAE-based modeling is best suited to the modeling of such systems: it enables a modeling based on first principles of the physics; it is physicsagnostic, and consequently accomodates arbitrary combinations of physics (mechanics, electrokinetics, hydraulics, thermodynamics, chemical reactions, etc.).

The pioneering work by Hilding Elmqvist [START_REF] Elmqvist | A structured model language for large continuous systems[END_REF] led to the emergence of the Modelica community in the 1990s, and the DAE-based modeling language of the same name (Modelica Association 2023) has become a de facto standard, with its object-oriented nature enabling a component-based modeling style. Its combined use with the port-Hamiltonian paradigm [START_REF] Rashad | Twenty years of distributed port-Hamiltonian systems: a literature review[END_REF] results in a methodology that is instrumental to the scalable modeling of large systems, additionally ensuring that the model architecture preserves the system architecture, in stark contrast to ODE-based modeling [START_REF] Benveniste | Multi-Mode DAE Models -Challenges, Theory and Implementation[END_REF]Benveniste, Caillaud, and Malandain 2022).

Consequently, DAE-based modeling requires that Modelica tools properly scale up to very large models. However, although Modelica enables the modeling of extremely large systems, its implementations (Dassault Systèmes 2002[START_REF]Modelica, A Unified Object-Oriented Language for Systems Modeling[END_REF][START_REF] Fritzson | The OpenModelica Integrated Environment for Modeling, Simulation, and Model-Based Development[END_REF] are often not capable of compiling and simulating such large models. Scaling has been and still is a subject for sustained effort by the Modelica community (Casella and Guironnet 2021), and although HPC issues belong to the landscape [START_REF] Braun | Solving large-scale Modelica models: new approaches and experimental results using OpenModelica[END_REF], a more specific issue is of uttermost importance for the Modelica language.

In the first steps of the compilation of a Modelica model, its hierarchical structure is flattened, thanks to a recursive syntactic inlining of the objects composing it. 1 The result of this flattening process is an unstructured DAE that can be exponentially larger than the source model. The structural analyses that are required for the generation of simulation code (namely, the index reduction of the DAE system, followed by a block-triangular form transformation of the reduced-index system) are then performed on this monolithic DAE model. As the compilation process does not fully take advantage of the hierarchical nature of the models it has to handle, the modeling capabilities offered by the Modelica language are undermined by performance issues on the structural analysis itself [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF][START_REF] Höger | Compiling Modelica: about the separate translation of models from Modelica to OCaml and its impact on variable-structure modeling[END_REF]). Additionally, model flattening poses a challenge when attempting to extend DAE-based modeling to higherorder modeling or dynamically changing systems [START_REF] Broman | Higher-Order Acausal Models[END_REF][START_REF] Broman | Meta-Languages and Semantics for Equation-Based Modeling and Simulation[END_REF][START_REF] Broman | Interactive Programmatic Modeling[END_REF].

In this paper, a new modular structural analysis algorithm is proposed that takes full advantage of the object tree structure of a DAE model. The bedrock of this method is a novel concept of structural analysis-aware interface for components. The essence of a component interface is to capture the necessary information about a Modelica class that needs to be exposed, in order to perform the structural analysis of a component comprizing instances of the former class, while hiding away useless information regarding the equations and all protected features it may contain.

In order to compute a component interface, one has to be able to perform the structural analysis of the possibly nonsquare DAE system that this component encapsulates, and to use the interfaces of the components it aggregates in this analysis. We base our algorithm on Pryce's Σ-method for index reduction [START_REF] Pryce | A Simple Structural Analysis Method for DAEs[END_REF], which essentially consists in the successive solving of two dual linear integer programs. The striking difference with Pryce's algorithm is that these problems are solved by parts, in a scalable manner.

Putting all of this together, it is then possible to perform a modular structural analysis, in which structural analysis is performed at the class level, and the results can then be instantiated for each component of the system model, knowing its context. Hence, structural information at the system level is derived from composing the result of component-level analysis. Modular structural analysis yields huge gains in terms of memory usage and computational costs, as the analysis of a single large-scale DAE is replaced with that of multiple smaller subsystems. Moreover, the analysis is performed at the class level, meaning that a single structural analysis is needed for all system components that are instances of the same class.

To the best of our knowledge, only [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF] addresses the specific issue of performing the structural analysis of a hierarchical model; Section 2.1 of this paper shows the approach proposed by Höger and explains why we regard it as an efficient but still partial solution. On the related subject of sorting equations, attention is paid in [START_REF] Zimmermann | Set-based graph methods for fast equation sorting in large DAE systems[END_REF] to methods that avoid unrolling loops and expanding arrays when sorting equations, an issue targeting the same overall objective as both [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF][START_REF] Höger | Compiling Modelica: about the separate translation of models from Modelica to OCaml and its impact on variable-structure modeling[END_REF]) and the present work.

Section 2 introduces two simple examples of DAE systems, to be used for illustrative purposes, and explains the modular approach from [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF] on one of them. In Section 3, we provide background information on structural analysis, which includes the Σ-method used to perform index reduction.

In Section 4, we introduce Σ-systems as an abstraction of DAE systems suitable for structural analysis. Σ-systems include the linear programming problems associated with the Σ-method and the structural description of the BTF. They also provide a notion of composition that abstracts the composition of DAE systems. This formalization leads us to the main contribution of this paper, presented in Section 5, where we propose the notion of Σ-interface. A Σ-interface exposes the necessary information to assemble partial structural analyses into a system-level structural analysis for a DAE system. We discuss how Σ-interfaces can reduce the computational cost of structural analysis for large systems.

Finally, Section 6 introduces our prototype implementation of the resulting modular method, then presents numerical applications to both examples presented in Section 2.

Examples

In this section, we develop two illustrative examples: a slightly modified version of the chained circuit proposed by C. Höger in [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF], and a homemade chained mass-and-spring system. Both will be used to illustrate how our approach changes and improves the structural analysis process of models where several instances of a same class are connected, while also highlighting different features of our algorithm.

A chained circuit

The circuit model

This example is a minor modification (with a second inductance added) of the one developed in [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF]. The following text is borrowed verbatim from this reference:

The Modelica representation of the circuit [...] consists of a name (Circuit), declarations of parameters (n), unknowns (u and i) and subcomponents (c). The physical behavior is defined directly by multiple equations, including the description of sub-circuit interconnection. Without knowledge about the internal structure of SubCircuit it is possible to use it inside the larger circuit, as long as it provides the corresponding interface (i.e. variables u and i). This composition of models is an easy and safe way to create more complex models out of simpler building blocks and is the very foundation of object-oriented modeling.

The Circuit and SubCircuit Modelica codes and schematics are given in Fig. 1 and 2 respectively.

The approach by C. Höger

While the flattening of the Circuit model is linear with n, the cost of the structural analysis of the resulting model is super-linear, thus preventing the classical approach from properly scaling up. In [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF], important contributions are proposed to cope with this problem. As far as we know, this paper is the first one pointing this issue very clearly. A faster method is proposed to perform the structural analysis of a hierarchical model, involving scoping and hiding.

This method is based on Pryce's Σ-method (Pryce 2001), which is explained in Section 3. The Σ-method involves solving a pair of dual Linear Programming problems (primal and dual), using a specific iterative algorithm for the dual. In Theorem 1 of [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF], the author provides a complexity argument to support the claim that the dual problem is not a bottleneck on its own.2 Hence, the author focuses on the primal problem, for which a decomposition method is proposed.

A recursive mass-spring-damper system

Figure 3 shows the model, consisting of a chain of massspring-damper elements, defined with the 1D translational components of the Modelica Standard Library (MSL). Although recursive classes are not allowed in Modelica, we use a recursive definition of a chain of elements. This gives a binary-tree structure to the model, that is best suited for the modular structural analysis method presented in the sequel of the paper. Our prototype implementation of the method supports recursive classes, with conditional statements evaluated statically, at compile time.

Our contribution for these examples

Despite its important contribution, we believe that [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF] does not provide the ultimate answer. While a hierarchical algorithm for solving the primal problem is a great contribution in terms of computational costs, it still does not fully take advantage of object-oriented modeling.

We would like instead to advance towards the separate compilation of components and systems, which consists in:

1. Proposing a notion of interface for model components, that is rich enough to subsequently perform systemlevel structural analysis; and 2. Proposing a modular structural analysis method, consisting of the needed algorithms performing hierarchi-cal structural analysis based on interface information. In doing so, both primal and dual problems, as well as the construction of a Block-Triangular Form (BTF) for the Jacobian, need to be addressed-this is also in contrast to the approach of [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF], where only the primal problem is considered.

To simplify our presentation, we skip the discussion of the BTF; it will be developed in an extended version of this paper, that is currently in the works.

The Σ-method for DAE

The index reduction method proposed by J. Pryce [START_REF] Pryce | A Simple Structural Analysis Method for DAEs[END_REF], called the Σ-method, is an interesting alternative to the classical method originally proposed by C. Pantelides [START_REF] Pantelides | The Consistent Initialization of Differential-Algebraic Systems[END_REF]. Its elegant and compact formulation as a pair of dual Linear Programming problems (LP) makes it particularly valuable for an extension to DAE components and architectures. The Σ-method can be summarized as follows:

The Σ-matrix For S = (F, X) a square DAE system involving equations f =0, where f ∈F, and variables x∈X and their derivatives, form the Σ-matrix Σ = (σ f,x ) ( f,x)∈F×X of the DAE system, where σ f,x is the highest differentiation order of variable x in f , or -∞ if x does not appear in f .

Primal problem

The primal LP encodes the search for a maximum weight transverse of Σ, that will be described as

( f , x f ) f ∈F or, equivalently, ( f x , x) x∈X in what follows.
The existence of a solution to the primal LP is a success check for the Σ-method.

Dual problem

The variables of the associated dual LP are variable offsets (d x ) x∈X and equation offsets (c f ) f ∈F , and we search for the minimal non-negative solution to this LP. [START_REF] Pryce | A Simple Structural Analysis Method for DAEs[END_REF] proves the uniqueness of this solution and proposes a relaxation method for finding it, given a solution to the primal problem.

The dual problem can be rewritten as the following constraint system, involving only the variable offsets

(d x ) x∈X : ∀x ∈ X : d x ≥ σ f x ,x ∀( f, x) ∈ E : d x -d x f ≥ σ f,x -σ f ,x f (1)
and the equation offsets are then given by

c f = d x f -σ f ,x f . (2)
It is proved in [START_REF] Pryce | A Simple Structural Analysis Method for DAEs[END_REF]) that the set of solutions of (1) does not depend on the particular choice of a solution to the primal problem.

Use of the offsets

Equation offset c f indicates how many times equation f =0 needs to be differentiated to get the index-reduced system, whereas d x indicates the maximum differentiation order of x in this index-reduced system. In addition, the solution of the primal problem is useful for computing a BTF for the Jacobian of this system. C. Höger [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF] explains that the primal problem is the main bottleneck, hence it focuses on solving it efficiently, by decomposing it into smaller subproblems. In contrast, we extend and adapt the Σ-method for open DAE systems, which are DAE systems that can possess more variables than equations, and that can be composed with other open DAE systems by unifying their common variables.

Structural analysis of open DAE

The set of variables of an open DAE system can be decomposed as X = X s ⊎X ℓ where X s and X ℓ are its sets of shared and local variables, respectively. For example, the SubCircuit of Fig. 2 possesses 2 shared variables and 10 local variables (declared as protected in Modelica).

Selectors

Since SubCircuit possesses 11 equations, we can regard it as a square system by assuming that one among the shared variables i, u is dependent and the other one is free (determined by the yet unspecified environment of this circuit):

X s = X free ⊎Y
where Y ⊆ X s is a subset of the shared variables. The set of dependent variables of the system is then

X dep = Y ⊎X ℓ .
If we compose open DAE system S with an environment S ′ (another open DAE system), then the two selectors Y and Y ′ for S and S ′ must satisfy the condition

Y ∩Y ′ = / 0 , (3)
expressing that S ′ cannot claim determining variable x if the latter is already claimed by S, i.e., belongs to Y . We say that selectors Y and Y ′ are compatible if (3) holds.

For SubCircuit, two possible choices for selectors are

Y = {i} or Y = {u} . (4)

Matching selectors in compositions

If {i} is selected, then u is free, i.e., it must be determined by the (yet unspecified) environment. Thus, we expect this environment to allow for a selector containing u but not i. This means that information (4) has to be exposed by S as part of its interface for structural analysis. Then, by exposing (4), S sets structural constraints on the interface of any environment for it. We will show that information (4) is sufficient for characterizing the environments that are compatible with S.

(5)

This information, however, is not sufficient to perform the structural analysis in a modular way. In the sequel, we will identify the information that is missing for this purpose.

Interfaces for the modular structural analysis of DAE systems

This section introduces the main contribution of this paper, which is the notion of component interface needed for the structural analysis of DAE systems. This notion is called Σ-interface as a reference to the Σ-method itself. Possible ways to perform the modular structural analysis of the Circuit example are also shown, as a way to illustrate the benefits of the modular approach.

The Σ-method for open DAE systems

We first show the result of the Σ-method, applied to the open DAE system SubCircuit (see Fig. 2), for the two possible choices of selectors given by (4).

Primal problem

The solutions of the primal problem for both selectors are given in Fig. 4. For each case, the free variable for the considered selector is written in blue. The chosen maximum weight transverse is indicated by highlighting in red the dependent variable associated to each equation, and we give in the third column the differentiation order of this variable in this equation. Two important observations will guide us in defining the primal Σ-interface:

1. The choice of the transverse, which fixes the assignment of variables to equations, does not depend on the (yet unspecified) environment of SubCircuit;

2. The total weight of the maximal transverse of SubCircuit, for a given selector, is added to the total weight of a maximal transverse of the environment (for a compatible selector), yielding the total weight of the overall transverse (that covers both the component and its environment).

Dual problem

The dual problem is the constraint system (1), whose dependent variables are the offsets (d x ) x∈X . For S = def SubCircuit, both selectors have to be considered:

• Note that, while edges ( f, x) are local (since all equations are local), variables can be shared between components; hence, dual problems only interact via the offsets of their shared variables (i and u for SubCircuit). This observation will guide us in defining the notion of dual Σ-interface.

Σ-interfaces

Primal Σ-interfaces

We say that a selector for a component is consistent if the primal problem has at least one solution for this selector.

Primal Σ-interface of a component Based on the discussion at the end of Section 5.1.1, it suffices to expose, for each consistent selector, the set of all pairs (selector, maximal transverse weight) (6) at the interface of the considered component, for the primal problem. Thus, (6) defines the primal Σ-interface of an open DAE system. For the SubCircuit example, the primal Σ-interface is the set ∀x, The primal Σ-interface of S is then obtained by collecting the pairs (Y, J Y S ) for every consistent selector of S. Given a consistent selector Y for S, we denote by

(Y = {i}, J Y ={i} S = 3) , (Y = {u}, J Y ={u} S = 2) (7)
d x ≥ 0 d uC , d iL1 , d iL2 ≥ 1 d iC -d uC ≥ -1 d uL1 -d iL1 ≥ -1 d uL2 -d iL2 ≥ -1 d u1 -d i1 ≥ 0 d i2 -d u2 ≥ 0 d u2 -d uL1 ≥ 0 d uC -d u1 ≥ 0 d u2 -d u1 ≥ 0 d i1 -d i2 ≥ 0 d iL1 -d i2 ≥ 0 d u1 -d uL2 ≥ 0 d uL1 -d uL2 ≥ 0 d u -d uL2 ≥ 0 d i1 -d iC ≥ 0 d iL2 -d iC ≥ 0 d iL2 -d i ≥ 0 ∀x, d x ≥ 0 d uC , d iL1 ≥ 1 d iC -d uC ≥ -1 d uL1 -d iL1 ≥ -1 d iL2 -d uL2 ≥ +1 d u1 -d i1 ≥ 0 d i2 -d u2 ≥ 0 d u2 -d uL1 ≥ 0 d uC -d u1 ≥ 0 d u2 -d u1 ≥ 0 d i1 -d i2 ≥ 0 d iL1 -d i2 ≥ 0 d u1 -d u ≥ 0 d uL1 -d u ≥ 0 d uL2 -d u ≥ 0 d i1 -d iC ≥ 0 d iL2 -d iC ≥ 0 d i -d iL2 ≥ 0
(π 1 (Y ), π 2 (Y ))
an optimizing selector pair. Remark that this pair may not be unique. However, the end result of the structural analysis (the offsets d x and c f ) does not depend upon the choice of (π 1 (Y ), π 2 (Y )).

Dual Σ-interfaces

Dual Σ-interface of a component Based on the discussion at the end of Section 5.1.2, we can now define the dual Σ-interface as collecting, for each consistent selector Y , the projection D Y of system (1) on the subset of offsets (d x ) x∈X s .

(9)

For the SubCircuit example, the dual Σ-interface is the following set (collecting two elements):

         Y = {i}, D Y ={i} : d i , d u ≥ 0 d i ≤ d u +1 , Y = {u}, D Y ={u} : d i , d u ≥ 0 d u ≤ d i +1          (10)
In ( 10), the two constraint systems are obtained by projecting, over the offsets d i and d u , the constraint systems given in Fig. 5.

Composing dual Σ-interfaces

Two open DAE systems S 1 and S 2 can only interact via their shared variables X s 1 ∪ X s 2 . Hence, for their composition S:

Given a consistent selector Y for S, the projection of the dual problem of S onto the offsets of

X s ⊆ X s 1 ∪ X s 2 is the composition of the projec- tions of the dual problems D π i (Y ) i of S i , i=1, 2 over the offsets of X s . (11)

Using Σ-interfaces in DAE systems

In this section, we illustrate the use of Σ-interfaces for the modular structural analysis of DAE systems. To each open DAE system, we associate its Σ-interface, by fusing the primal and dual Σ-interfaces defined above: it is a set of triples whose elements are a consistent selector, the weight of a solution of the corresponding primal problem, and the projection of the dual problem on the offsets of the shared variables. We shall then use ( 8) and ( 11) to perform the structural analysis of a DAE system in a modular way.

The Circuit, seen as a chain

From ( 7) and ( 10), one gets the following Σ-interface for SubCircuit:

         Y = {i} , J Y S = 3, d i , d u ≥ 0 d i ≤ d u +1 Y = {u} , J Y S = 2, d i , d u ≥ 0 d u ≤ d i +1          (12)
In what follows, successive instances of SubCircuit are chained, as shown in Fig. 1. The structural analysis will be performed by induction on the length n of the chain.

Call S n the chain of length n; in particular, S 1 = SubCircuit, denoted by SC in what follows. Call i, u n the shared variables of S n . Seen an a chain, S n is obtained by composing S n-1 with SC and adding a Kirchhoff equation for voltages. We regard this last equation as a component with no local variables, denoted by eq below.

Note that we also have to rename u as v in SC in order to avoid name clashes, which we write [u/v], and that variable hiding has to be used on the result of the composition, as variables u n-1 and v have to be made local in S n .

As a result, for all n ≥ 2, S n is defined as the following composition:

S n = hide u n-1 , v in    S n-1 SC[u/v] u n = u n-1 + v (13)
For every n, S n has exactly one more variable than it has equations. Thus, there are two possible selectors for S n : Y S = {i} and Y S n = {u n }.

Case Y S n = {u n } Then, i is free and there is only one triple of compatible selectors in the composition occurring in the right-hand side of (13):

Y S n-1 = {u n-1 },Y SC = {v},Y eq = {u n } .
An immediate induction argument shows that the optimal weight for S n is 2(n-1) + 0 + 2 = 2n.

We prove by induction that the dual Σ-interface is

d i , d u n ≥ 0 d u n ≤ d i + 1 (14)
for every n. This holds for n = 1, as ( 14) then yields the dual Σ-interface of class SubCircuit for selector {u} (see ( 10)). Assuming that ( 14) holds for n-1, the dual Σ-interfaces compose as follows:

hide u n-1 , v in            d i , d u n-1 , d v , d u n ≥ 0 d u n-1 ≤ d i + 1 d v ≤ d i + 1 d u n ≤ d u n-1 d u n ≤ d v
which yields ( 14) when projected on u n , i.

Consequently, the variable offsets of a SubCircuit in a Circuit are independent of the number of chained instances, and the equation offset of an equation in the k-th component of the chain is equal to that of the same equation in the original class.

Case Y S n = {i} Then, u n is free and there are two triples of compatible selectors:

Y S n-1 = {u n-1 },Y SC = {i},Y eq = {v} (15) Y S n-1 = {i},Y SC = {v},Y eq = {u n-1 } (16)
By adding the contributions of the components, in the order they are written above, in each case, we get:

• In case ( 15): [2(n-1)] + 0 + 3 = 2n + 1 • In case ( 16): [2(n-1) + 1] + 0 + 2 = 2n + 1 Thus, each of the two triples ( 15) and ( 16) is optimizing, which brings an important question: Does the dual Σ-interface depend on the choice of a tuple of selectors? If triple ( 15) is used: System (14) for n-1 provides us with the dual Σ-interface of S n-1 for selector u n-1 . For S n , the dual Σ-interfaces then compose as follows:

hide u n-1 , v in            d i , d u n-1 , d v , d u n ≥ 0 d u n-1 ≤ d i + 1 d i ≤ d v + 1 d v ≤ d u n-1 d v ≤ d u n which yields d i , d u n ≥ 0 d i ≤ d u n + 1 (17)
If triple ( 16) is used: We will prove by induction that (17) still holds. Assuming that (17) holds for n-1, the dual Σ-interfaces compose as follows:

hide u n-1 , v in            d i , d u n-1 , d v , d u n ≥ 0 d i ≤ d u n-1 + 1 d v ≤ d i + 1 d u n-1 ≤ d v d u n-1 ≤ d u n
which yields again (17).

Recall that solutions to the Σ-method's dual problem are independent of the particular choice of a solution for the primal problem. This property generalizes to the composition of dual interfaces: the dual Σ-interface of a composition does not depend on the choice of a particular optimizing tuple.

(18)

The Circuit, seen as a tree

Instead of a chain, the architecture of the Circuit can be regarded as a binary tree:

T m = def ŜC ↙ ↘ T - m-1 T + m-1 (19)
where:

• the length ℓ m of T m is defined by

ℓ m = 2ℓ m-1 + 1, which is exponential in m; • T m has shared variables i m , u m ; • T ± m-1 are two copies of T m-1 with renamings [i ± m-1 /i];
• ŜC is made of an instance of SC and a component (similar to component eq above) made of the connection equation

u m = u - m-1 + u + u + m-1 .
We reuse the Σ-interface (12) of class SubCircuit, the optimal weights 2ℓ m and 2ℓ m +1 for selectors Y T m ={u m } and Y T m ={i} respectively, and we focus on the dual interface. It is computed by the recursion shown in Fig. 6.

Y T m ={u m } : hide u, u ± m-1 in                      0 ≤ d i , d u ± m-1 , d u m d u - m-1 ≤ d i + 1 d u + m-1 ≤ d i + 1 d u ≤ d i + 1 d u m ≤ d u ± m-1 d u m ≤ d u which yields 0 ≤ d i , d u m d u m ≤ d i + 1 Y T m ={i} : hide u, u ± m-1 in                      0 ≤ d i , d u ± m-1 , d u m d u - m-1 ≤ d i + 1 d u + m-1 ≤ d i + 1 d i ≤ d u + 1 d u ≤ d u ± m-1 d u ≤ d u m which yields 0 ≤ d i , d u m d i ≤ d u m + 1 Figure 6
. The dual Σ-interface of tree shaped architecture (19).

Discussion

It is worth comparing the chain-based and tree-based approaches above. At first glance, since recursion arguments were used in both cases, the two approaches may seem equivalent in terms of computational costs. However, this impression shall not last once we detail how implementations should proceed.

• For the chain architecture of Section 5.3.1, each induction step consists in the computation of the Σ-interface of prefix S k as a function of the Σ-interface of S k-1 , for k increasing from 2 to n.

• For the tree-shaped architecture of Section 5.3.2, the induction step expresses the Σ-interface of the root of each subtree T k as a function of the Σ-interface of each subtree T ± k-1 . Remark that the Σ-interface of each subtree T k is only computed once, as all leaf components are instances of the same class, and all subtrees of height k have, by construction, identical interfaces.

The number of steps is proportional to n in the first case, and m ∼ log n in the second case. Since the computational complexity of each induction step is roughly the same, treating the Circuit model as a tree-shaped architecture can dramatically improve performance.

6 Implementation and experimental results

Implementation considerations

A significant difficulty in our modular structural analysis is the consideration of selector-dependent structural analyses.

Both the primal and dual problems are functions of the selectors, which can be numerous for components with a large number of public variables. Although Modelica models are usually sparse, with components exposing only a few variables, a mere enumeration of selectors can result in an exponential growth of the handled data structures.

To deal with this issue, we advocate the approach proposed in [START_REF] Benveniste | Algorithms for the Structural Analysis of Multimode Modelica Models[END_REF]Caillaud, Malandain, and Thibault 2020) for multimode DAE systems. In these works, a dual representation is introduced, where equations are labeled with a predicate on mode variables, characterizing the modes under which they are active. Using this representation (instead of a direct one, that lists the active equations for each mode of the model) provides a compact representation of the structure of multimode systems. Reduced Ordered Binary Decision Diagrams, or ROBDD [START_REF] Bryant | Graph-Based Algorithms for Boolean Function Manipulation[END_REF], provide not only an adapted data structure, but also efficient computation algorithms that can be used for performing the whole structural analysis in an "all-modes-at-once" fashion.

For modular structural analysis, a similar representation can be used, but with selectors instead of modes. The whole structural analysis chain can then be performed in an "all-selectors-at-once" fashion: [START_REF] Benveniste | Algorithms for the Structural Analysis of Multimode Modelica Models[END_REF] provides all building blocks for the computation of primal interfaces and their compositions, as presented in Section 5.2.1. The computation and composition of dual interfaces reduce to projecting parametrized constraint systems such as the one illustrated Fig. 5. Such systems are called Difference Bound Matrices (DBM) and come equipped with a rich calculus [START_REF] Dill | Timing Assumptions and Verification of Finite-State Concurrent Systems[END_REF][START_REF] Miné | A New Numerical Abstract Domain Based on Difference-Bound Matrices[END_REF]) with a polynomial computational complexity. Although there are excellent implementations of DBM, parametric DBM remain to be implemented, possibly using tuples of ROBDD for the representation of matrix elements.

Benchmarks and measured performance

Experimental results have been obtained on the modular structural analysis method, using two benchmarks: the chain circuit (Section 2.1, Figure 2) and the mass-springdamper model (Section 2.2, Figure 3). Models of increasing sizes have be analyzed, up to n = 10 12 . The experimental results presented below have been obtained with a prototype implementation of the method, based on an enumeration of the selectors. This has a limited impact on performance, for these particular models, since the selector combinatorics is limited to 2 cases for the chain circuit, and 6 cases for the mass-spring-damper model. The key feature of the prototype implementation is the use of a dynamic programming approach for the computation of the interfaces of model tree nodes. For this purpose, a memoization table (see [START_REF] Cormen | Introduction to algorithms[END_REF], pages 390-392) is used to store the computed interfaces. The software consists in about 10 kLOC of OCaml code and performance has been measured on a MacBook Pro with a 2.4GHz 8core i-9 Intel processor with 16GB of RAM.

Figure 7 shows the computation times for both models as a function of parameter n. It clearly appears that the empirical time complexity of the method is a logarithmic function of n, like in [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF]. Memory usage is very modest, with about 1MB to store the memoization table.

Conclusion and perspectives

In this paper, we present a modular index-reduction method for (possibly incomplete) DAE systems, based on John Pryce's Σ-method [START_REF] Pryce | A Simple Structural Analysis Method for DAEs[END_REF]) and extending the seminal work of Christoph Höger on scalable algorithms for the compilation of Modelica [START_REF] Höger | Faster Structural Analysis of Differential-Algebraic Equations by Graph Compression[END_REF][START_REF] Höger | Compiling Modelica: about the separate translation of models from Modelica to OCaml and its impact on variable-structure modeling[END_REF]. Our method is built upon three key contributions: (i) a concept of selector that allows to characterize, from a structural analysis point of view, the possible effects of unknown environments on an incomplete DAE system; (ii) a concept of interface for the primal and dual problems of the Σ-method, that encapsulates the minimal information regarding a subsystem that needs to be exposed to its environment in order to perform the structural analysis; and (iii) interface composition and transformation operators that allow to compute the interface of a system from the interfaces of its parts.

Our modular structural analysis method is well-suited to the Modelica language, since the interface of a class can be computed inductively from the interfaces of the objects contained in the class. Modelica models are often sparse, meaning that each component shares only a few variables with its environment. This guarantees that the interface of a component remains small, independently of the number of components, variables and equations it may contain.

We believe that this concept of interface, and the modular structural analysis method it yields, pave the way towards a genuine separate compilation of Modelica, that can scale up to extremely large models.

The benchmarks performed with our prototype implementation demonstrate that extremely large models, organized in a component tree with sufficient regularity, can be analyzed in a few milliseconds.

Future work shall focus on a more robust implementation of the method, based on our IsamDAE multimode DAE structural analysis software [START_REF] Benveniste | Algorithms for the Structural Analysis of Multimode Modelica Models[END_REF]Caillaud, Malandain, and Thibault 2020). We plan to use a functional, BDD-based, representation of interfaces, in order to curb the combinatorics of selectors that is expected when computing interfaces of Modelica classes with a large number of public variables. This functional representation is also a promising approach to the extension of the notion of interface to multimode models.

The structural analysis of Modelica models comprising for loops could be done by (i) computing a tree decomposition of the component graph obtained from the evaluation of the loops, followed by (ii) the modular structural analysis of the resulting component tree. Benchmarking this approach on the scalable test suite (Casella and Guironnet 2021) would be of great interest.

Several features of the Modelica language, such as inner/outer declarations, expandable connectors and overconstrained connections, may turn out to be difficult to deal with. In the Modelica Language Specification (Modelica Association 2023), the semantics of these features is expressed in terms of an elaboration phase, that is not modular. Devising a modular transformation of these features into the Modelica kernel language is a challenge that needs to be addressed before the modular structural analysis method can be applied to the full Modelica language.

There are also a few fundamental issues, related to the modular structural analysis, that ought to be investigated: (i) Is it possible to decide, by inductive reasoning, whether a parametric model is structurally nonsingular for every valuation of its parameters? (ii) Is it possible to perform the structural analysis of DAE systems that have a large treewidth, such as systems organized as a grid of dimension 2 or higher?

Answers to these questions are key to the design of scalable separate compilation methods for the Modelica language. One could envision the standardization of interfaces, possibly as an extension of the FMI standard for model exchange, so that precompiled Modelica libraries, equipped with their Σ-interfaces, could be built and reused.

Figure 2 .

 2 Figure 2. Modelica code of the SubCircuit class, and schematics of the corresponding electrical circuit.

Figure 3 .

 3 Figure 3. A mass-spring-damper system; the element (top-left) is assembled from the 1-D translational mechanical components of the Modelica Standard Library. An assembly of two elements is shown at the bottom-left. A chain of mass-spring-damper elements of length n is defined by the recursive Modelica-like class shown on the right. Although Modelica does not allow for recursive classes, our software prototype allows recursion, provided conditional statements can be evaluated at compile-time.

  S has selector {i}. A maximal weight transverse is shown in red on Fig. 4-left. With this transverse, the dual problem (1) is shown in Fig. 5-left. • S has selector {u}. A maximal weight transverse is shown in red on Fig. 4-right. With this transverse, the dual problem (1) is shown in Fig. 5-right.

Figure 4 .

 4 Figure 4. The primal problem for SubCircuit, seen as an open DAE system, for two possible choices for the selector. Left: selector Y = {i} yields a maximal transverse (in red) of weight 3; right: selector Y = {u} yields a maximal transverse (in red) of weight 2. In each case, the contribution of each equation to the weight of the transverse is provided on the right of the equation, and the free variable (to be determined by the environment) is highlighted in blue.

Figure 5 .

 5 Figure 5. The dual problem when S has selector {i} (left) and {u} (right).

Figure 7 .

 7 Figure 7. Number of interfaces computed (top plot) and computation time (bottom plot, in ms) for the mass-spring-damper (purple curves) and chain circuit (green curves) models, as functions of parameter n.

model Circuit parameter Integer n = 10; SubCircuit[n] c; Real u,i; equation c[1].i = i; u = sin(time); for j in 2:n loop c[j].i = c[j-1].i; end for; sum(c.u) = u; end Circuit;

	c_1	c_2	...	c_n
		u
	Figure 1. Modelica code and electrical schematics of the chained circuit.

class SubCircuit parameter Real R1, R2, L, C; Real i, u; protected Real u1,i1, u2,i2, uC,iC, uL1,iL1, uL2,iL2; equation iC=C * der(uC); uL1=L1 * der(iL1); uL2=L2 * der(iL2); u1 = R1 * i1; u2 = R2 * i2; u2 = uL1; uC=u1+u2; i1=i2+iL1; u=u1+uL1+uL2; i=i1+iC; i=iL2; end SubCircuit;

See (Modelica Association

2023), Section 5.6 for a complete definition of the flattening process.

Our own numerical experiments confirm this observation.