Roland Benoit
email: roland.benoit@cnrs-orleans.fr

Nesrine Belhadj
email: nesrine.belhadj@cnrs-orleans.fr

Maxence Lailliau
email: maxence.lailliau@cnrs-orleans.fr

Philippe Dagaut
email: dagaut@cnrs-orleans.fr

Autoxidation of terpenes, a common pathway in tropospheric and low temperature combustion conditions: the case of limonene and α-pinene

Introduction

The links between atmospheric and combustion chemistry have often been studied from the point of view of tropospheric reactions of combustion effluents or pollutants, e.g., oxidation of volatile organic compounds, nitrogen oxides reactions, sulfur chemistry [START_REF] Barsanti | Formation of Low-Volatility Organic Compounds in the Atmosphere: Recent Advancements and Insights[END_REF][START_REF] Shrivastava | Recent advances in understanding secondary organic aerosol: Implications for global climate forcing[END_REF][START_REF] Kenseth | Synergistic O3 + OH oxidation pathway to extremely low-volatility dimers revealed in β-pinene secondary organic aerosol[END_REF][START_REF] Bianchi | Highly Oxygenated Organic Molecules (HOM) from Gas-Phase Autoxidation Involving Peroxy Radicals: A Key Contributor to Atmospheric Aerosol[END_REF]. Climate change and the increase of large wildfire events have profoundly modified the relationship between atmospheric chemistry and combustion at large. When studying the oxidation of chemicals and the formation of SOA in the atmosphere, it becomes necessary to determine the contribution of different oxidation pathways pertaining to atmospheric chemistry, combustion chemistry, or both.

In low-temperature combustion (cool flame) the formation of oxidized organic molecules (OOM) is mainly attributed to autoxidation reactions [START_REF] Belhadj | Experimental and kinetic modeling study of n-hexane oxidation. Detection of complex low-temperature products using high-resolution mass spectrometry[END_REF][START_REF] Benoit | On the similarities and differences between the products of oxidation of hydrocarbons under simulated atmospheric conditions and cool flames[END_REF], whereas in atmospheric chemistry, it is only relatively recently that this pathway has been considered [START_REF] Vereecken | Low-volatility poly-oxygenates in the OH-initiated atmospheric oxidation of α-pinene: impact of non-traditional peroxyl radical chemistry[END_REF][START_REF] Crounse | Autoxidation of organic compounds in the atmosphere[END_REF]Jokinen et al., 2014a;[START_REF] Berndt | Gas-Phase Ozonolysis of Cycloalkenes: Formation of Highly Oxidized RO2 Radicals and Their Reactions with NO, NO2, SO2, and Other RO2 Radicals[END_REF]Jokinen et al., 2015a;[START_REF] Berndt | Hydroxyl radical-induced formation of highly oxidized organic compounds[END_REF][START_REF] Iyer | Molecular mechanism for rapid autoxidation in α-pinene ozonolysis[END_REF]. Also, it has been identified that highly oxygenated molecules (HOMs), a source of secondary organic aerosols (SOA), can result from autoxidation processes [START_REF] Wang | Efficient alkane oxidation under combustion engine and atmospheric conditions[END_REF][START_REF] Tomaz | Structures and reactivity of peroxy radicals and dimeric products revealed by online tandem mass spectrometry[END_REF][START_REF] Bianchi | Highly Oxygenated Organic Molecules (HOM) from Gas-Phase Autoxidation Involving Peroxy Radicals: A Key Contributor to Atmospheric Aerosol[END_REF]. Modeling studies complemented by laboratory experiments showed that autoxidation mechanisms proceed simultaneously on different RO2˙ radicals and leads to, through isomerization and addition of O2, the production of a wide range of oxidized compounds in a few hundredths of a second (Jokinen et al., 2014a;[START_REF] Berndt | Hydroxyl radical-induced formation of highly oxidized organic compounds[END_REF][START_REF] Bianchi | Highly Oxygenated Organic Molecules (HOM) from Gas-Phase Autoxidation Involving Peroxy Radicals: A Key Contributor to Atmospheric Aerosol[END_REF]. Autoxidation is based on an H-shift and oxygen addition which starts with the initial production of RO2˙ radicals. This mechanism can repeat itself several times and lead to recurrent oxygen additions to form HOMs: R+ O2 ⇄ ROO • (first O2-addition); ROO • ⇄ • QOOH (H-shift); • QOOH +O2 ⇄ • OOQOOH (second O2-addition);

• OOQOOH ⇄ HOOQ' • OOH (H-shift); HOOQ' • OOH +O2 ⇄ (HOO)2Q'OO • (third O2-addition); (HOO)2Q'OO • ⇄ (HOO)2Q" • OOH (H-shift); (HOO)2Q" • OOH +O2 ⇄ (HOO)3Q" OO • (fourth O2-addition) etc. Recent works have shown that, under certain atmospheric conditions, this autoxidation mechanism could be competitive with other reaction pathways involving RO2˙ radicals [START_REF] Bianchi | Highly Oxygenated Organic Molecules (HOM) from Gas-Phase Autoxidation Involving Peroxy Radicals: A Key Contributor to Atmospheric Aerosol[END_REF], e.g., the carbonyl channel (ROO • → R-HO + OH), the hydroperoxide channel (ROO • + HOO • → ROOH + O2 and RO • + • OH + O2), disproportionation reactions (ROO • + R'OO • → RO • + R'O • + O2 and R-HO + R'OH + O2), accretion reactions (ROO • + R'OO • → ROOR' + O2). A study published in 2021 showed that the oxidation of alkanes follows this autoxidation mechanism under both atmospheric and combustion conditions [START_REF] Wang | Efficient alkane oxidation under combustion engine and atmospheric conditions[END_REF]. Also, that work confirmed that internal H-shifts in autoxidation can be promoted by the presence of functional groups, as predicted earlier [START_REF] Otkjaer | Calculated Hydrogen Shift Rate Constants in Substituted Alkyl Peroxy Radicals[END_REF] for RO2˙ radicals containing OOH, OH, OCH3, CH3, C=O, or C=C groups. To further assess the importance of these pathways, available data must be compared along with their experimental specificities. In laboratory studies under simulated atmospheric conditions, oxidation occurs at near-ambient temperatures (250-300 K), at atmospheric pressure, in the presence of ozone and/or OH˙ radicals, and with low initial terpene concentrations. In combustion, the OH˙ radical, temperature, and pressure are driving autoxidation. Initial reactant concentrations are generally higher compared to atmospheric conditions, so as to compensate for the absence of ozone and get oxidation to proceed, since terpenes, as other hydrocarbons, react very slowly with O2. Rising temperature increases isomerization rates and favors autoxidation, at the expense of other possible reactions of RO2˙ radicals. In atmospheric chemistry, at near room temperature, autoxidation can be initiated via two reaction mechanisms: reactions with ozone or with the OH˙ radical. Suppression of one of these pathways by scavengers https://doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.

generally changes the amount of products formed (sometimes a drastic decrease >90% was observed [START_REF] Kenseth | Synergistic O3 + OH oxidation pathway to extremely low-volatility dimers revealed in β-pinene secondary organic aerosol[END_REF][START_REF] Meusinger | Chemical and isotopic composition of secondary organic aerosol generated by α-pinene ozonolysis[END_REF][START_REF] Pospisilova | On the fate of oxygenated organic molecules in atmospheric aerosol particles[END_REF]), but does not affect the diversity of observed chemical formulae [START_REF] Kenseth | Synergistic O3 + OH oxidation pathway to extremely low-volatility dimers revealed in β-pinene secondary organic aerosol[END_REF][START_REF] Pospisilova | On the fate of oxygenated organic molecules in atmospheric aerosol particles[END_REF]. It has been reported earlier that a temperature rise from 250 to 273K does not affect the distribution of HOMs [START_REF] Quéléver | Effect of temperature on the formation of highly oxygenated organic molecules (HOMs) from alpha-pinene ozonolysis[END_REF] whereas Trostl et al. suggested that the distribution of HOMs is affected by temperature, α-pinene or particle concentration [START_REF] Tröstl | The role of low-volatility organic compounds in initial particle growth in the atmosphere[END_REF]. Similarly, the experiments of Huang et al. performed at different temperatures (223 K and 296 K) and precursor concentration (α-pinene 0.714 and 2.2 ppm) suggested that the physicochemical properties, such as the composition of the oligomers, can be affected by a variation of temperature [START_REF] Huang | α-Pinene secondary organic aerosol at low temperature: chemical composition and implications for particle viscosity[END_REF]. Both in combustion and atmospheric chemistry, autoxidation can yield highly oxygenated molecules (HOMs), e.g., compounds containing more than 7 oxygen atoms [START_REF] Benoit | On the similarities and differences between the products of oxidation of hydrocarbons under simulated atmospheric conditions and cool flames[END_REF][START_REF] Bianchi | Highly Oxygenated Organic Molecules (HOM) from Gas-Phase Autoxidation Involving Peroxy Radicals: A Key Contributor to Atmospheric Aerosol[END_REF]. The term 'HOM' is generally associated with atmospheric chemistry [START_REF] Bianchi | Highly Oxygenated Organic Molecules (HOM) from Gas-Phase Autoxidation Involving Peroxy Radicals: A Key Contributor to Atmospheric Aerosol[END_REF], but this nomenclature does not specify the chemical properties of a compound.

In other words, in combustion we can also observe highly oxidized chemical compounds similar to those relevant to atmospheric chemistry. The broad range of chemical molecules formed and the impact of the experimental conditions on their character remain a subject for atmospheric chemistry as well as combustion chemistry studies.

Moreover, whatever the mechanism of aerosol formation, i.e., oligomerization, addition, or accretion, their composition will be linked to that of the initial radical pool [START_REF] Tomaz | Structures and reactivity of peroxy radicals and dimeric products revealed by online tandem mass spectrometry[END_REF].

In low-temperature combustion, when the temperature is increased, autoxidation rate goes through a maximum between 500 and 670 K, depending on the nature of the fuel [START_REF] Belhadj | Oxidation of di-n-butyl ether: Experimental characterization of low-temperature products in JSR and RCM[END_REF][START_REF] Belhadj | Experimental and kinetic modeling study of n-hexane oxidation. Detection of complex low-temperature products using high-resolution mass spectrometry[END_REF], and the rate of formation of HOM also increases [START_REF] Bianchi | Highly Oxygenated Organic Molecules (HOM) from Gas-Phase Autoxidation Involving Peroxy Radicals: A Key Contributor to Atmospheric Aerosol[END_REF]. Therefore, temperature has a measurable effect on the amount of HOMs formed [START_REF] Wang | Efficient alkane oxidation under combustion engine and atmospheric conditions[END_REF]. In low-temperature combustion chemistry as in atmospheric chemistry, the oxidation of a chemical compound leads to the formation of several thousands of chemical products which result from successive additions of oxygen, isomerization, accretion, fragmentation, and oligomerization [START_REF] Benoit | On the similarities and differences between the products of oxidation of hydrocarbons under simulated atmospheric conditions and cool flames[END_REF]. The exhaustive analysis of chemical species remains, in the current state of instrumental limitations, impossible. This would consist in analyzing several thousands of molecules in MS² mode using separative techniques such as ultra-high pressure liquid chromatography (UHPLC) or ion mobility spectrometry (IMS) [START_REF] Krechmer | Ion mobility spectrometry-mass spectrometry (IMS-MS) for on-and offline analysis of atmospheric gas and aerosol species[END_REF][START_REF] Kristensen | High-Molecular Weight Dimer Esters Are Major Products in Aerosols from α-Pinene Ozonolysis and the Boreal Forest[END_REF]. Nevertheless, it is possible to classify these molecular species, considering only CxHyOz compounds, according to various criteria accessible via graphic tools representation such as van Krevelen diagrams, double bond equivalent number (DBE), and average carbon oxidation state (OSc) versus the number of carbon atoms [START_REF] Kourtchev | Molecular Composition of Fresh and Aged Secondary Organic Aerosol from a Mixture of Biogenic Volatile Compounds: A High-Resolution Mass Spectrometry Study[END_REF][START_REF] Nozière | The Molecular Identification of Organic Compounds in the Atmosphere: State of the Art and Challenges[END_REF]. Such postprocessing of large datasets has the advantage of immediately highlighting families of compounds or physicochemical properties such as the condensation of molecules (vapor pressure), the large variety of oxygenated products (CxHyO1-15 in the present experiments) and the formation of oligomers [START_REF] Kroll | Carbon oxidation state as a metric for describing the chemistry of atmospheric organic aerosol[END_REF]Xie et al., 2020).

In addition to the recent studies focusing on the first steps of autoxidation, a more global approach, based on the comparison of possible chemical transformations related to autoxidation in low temperature combustion and atmospheric chemistry, is needed for evaluating the importance of autoxidation under tropospheric and lowtemperature combustion conditions. In order to study the effects of ozonolysis and combustion on the diversity of chemical molecules formed by autoxidation, we have selected α-pinene and limonene, two terpene isomers among the most abundant in the troposphere. Limonene has a single ring structure and two double bonds, one of which is exocyclic. α-Pinene has a bicyclic structure and a single endo-cyclic double bond. These two isomers with their distinctive characters are good candidates for studying autoxidation versus initial chemical structure and https://doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License. temperature. For α-pinene, in addition to the reactivity of its endo-cyclic double bond, products of ring opening of the cyclobutyl group have been detected, which could explain the diversity of observed oxidation products [START_REF] Kurtén | Computational Study of Hydrogen Shifts and Ring-Opening Mechanisms in α-Pinene Ozonolysis Products[END_REF][START_REF] Iyer | Molecular mechanism for rapid autoxidation in α-pinene ozonolysis[END_REF]. This large pool of products is increased in the case of limonene by the presence of two double bonds [START_REF] Hammes | Carboxylic acids from limonene oxidation by ozone and hydroxyl radicals: insights into mechanisms derived using a FIGAERO-CIMS[END_REF]Jokinen et al., 2015a). Nevertheless, the diversity of oxidation products can be compensated by similarities in terms of reaction mechanisms and intermediates formed through autoxidation [START_REF] Savee | Carbon radicals. Direct observation and kinetics of a hydroperoxyalkyl radical (QOOH)[END_REF].

In the present work, which is a prolongation of that published earlier for the oxidation of limonene alone [START_REF] Benoit | On the similarities and differences between the products of oxidation of hydrocarbons under simulated atmospheric conditions and cool flames[END_REF], we oxidized α-pinene and limonene in a jet-stirred reactor at atmospheric pressure, excess of oxygen, and elevated temperature. Then, we characterized the impact of using different ionization techniques (HESI and APCI) in positive and negative modes on the pool of detected chemical formulae. The particularities of each ionization mode were analyzed to identify the most suitable ionization technique for exploring the formation of autoxidation products under low temperature combustion. Chemical formulae detected here and in atmospheric chemistry studies were compiled and tentatively used to evaluate the importance of autoxidation under both conditions.

Experiments

Oxidation experiments

The present experiments were carried out in a fused silica jet-stirred reactor (JSR) setup presented earlier [START_REF] Dagaut | A jet-stirred reactor for kinetic studies of homogeneous gas-phase reactions at pressures up to ten atmospheres (≈1 MPa)[END_REF]Dagaut et al., 1988) and used in previous studies [START_REF] Dagaut | Kinetic Modeling of Propane Oxidation[END_REF][START_REF] Benoit | On the similarities and differences between the products of oxidation of hydrocarbons under simulated atmospheric conditions and cool flames[END_REF][START_REF] Belhadj | Experimental and kinetic modeling study of n-hexane oxidation. Detection of complex low-temperature products using high-resolution mass spectrometry[END_REF]. The two isomers α-pinene and limonene were studied separately. As in earlier works [START_REF] Benoit | On the similarities and differences between the products of oxidation of hydrocarbons under simulated atmospheric conditions and cool flames[END_REF][START_REF] Belhadj | Experimental and kinetic modeling study of n-hexane oxidation. Detection of complex low-temperature products using high-resolution mass spectrometry[END_REF], α-pinene (+), 98% pure from Aldrich and limonene (R)-(+), >97% pure from Sigma, were pumped by an HPLC pump (Shimadzu LC10 AD VP) with an online degasser (Shimadzu DGU-20 A3) and sent to a vaporizer assembly where it was diluted by a nitrogen flow. Each terpene isomer and oxygen, both diluted by N2, were sent separately to a 42 mL JSR to avoid oxidation before reaching 4 injectors (nozzles of 1 mm I.D.)

providing stirring. The flow rates of nitrogen and oxygen were controlled by mass flow meters. Good thermal homogeneity along the vertical axis of the JSR was recorded (gradients of < 1 K/cm) by thermocouple measurements (0.1 mm Pt-Pt/Rh-10% wires located inside a thin-wall silica tube). In order to observe the oxidation of these isomers, which are not prone to strong self-ignition, the oxidation of 1% of these chemical compounds (C10H16) under lean fuel conditions (equivalence ratio 0.25, 56% O2, 43% N2), experiments were carried out at 590 K, atmospheric pressure, and a residence time of 2 s. Under these conditions, the oxidation of the two isomers is initiated by slow H-atom abstraction by molecular oxygen. The fuel radicals react rapidly with O2 to form peroxy radicals which undergo further oxidation, characteristic of autoxidation. The absence of ozone, and no need for the addition of a scavenger, allow probing reaction mechanisms and observing chemical species potentially specific to autoxidation initiated by OH radicals. A 2 mm I.D. probe was used to collect samples. To measure lowtemperature oxidation products ranging from early oxidation steps to highly oxidized molecules, the sonic probe samples were bubbled in cooled acetonitrile (UHPLC grade ≥99.9, T= 0°C, 250 mL) for 90 min. The resulting solution was stored in a freezer at -15°C.

Chemical analyses

Analyses of samples collected in acetonitrile were performed by direct sample instillation (rate: 3μL/min and recorded for 1 min for data averaging) in the ionization chamber of a high-resolution mass spectrometer (Thermo Scientific Orbitrap® Q-Exactive, mass resolution 140,000 and mass accuracy <0.5 ppm RMS). Both heated electrospray ionization (HESI) and atmospheric chemical ionization (APCI) were used in positive and negative modes for the ionization of products. HESI settings were: spray voltage 3.8 kV, T vaporizer of 150°C, T capillary 200°C, sheath gas flow of 8 arbitrary units (a.u.), auxiliary gas flow of 1 a.u., sweep gas flow of 0 a.u.. In APCI, settings were: spray voltage 3.8 kV, vaporizer temperature of 150°C, capillary temperature of 200°C, sheath gas flow of 8 a.u., auxiliary gas flow of 1 a.u., sweep gas flow of 0 a.u., corona current of 3μA. In order to avoid transmission and detection effects of ions depending on their mass inside the C-Trap [START_REF] Hecht | Fundamentals and Advances of Orbitrap Mass Spectrometry[END_REF], acquisitions with three mass ranges were performed (m/z 50-750; m/z 150-750; m/z 300-750). The upper limit of m/z 750 was chosen because of the absence of signal beyond this value. We verified that no significant oxidation occurred in the HESI and APCI ion sources by injecting a limonene-ACN mixture. The optimization of the Orbitrap ionization parameters in HESI and APCI did not show any clustering phenomenon for these two isomers.

The parameters evaluated were: injection source -capillary distance, vaporization and capillary temperatures, applied difference of potential, injected volume, flow rate of nitrogen in the ionization source. Positive and negative HESI mass calibrations were performed using Pierce TM calibration mixtures (Thermo Scientific).

Chemical compounds with relative intensity less than 1 ppm to the highest mass peak in the mass spectrum were not considered. Nevertheless, it should be considered that some of the molecules presented in this study could result from our experimental conditions (continuous flow reactor, reagent concentration, temperature, reaction time) and to some extent from our acquisition conditions, different from those of the previous studies [START_REF] Deng | Temperature and acidity dependence of secondary organic aerosol formation from α-pinene ozonolysis with a compact chamber system[END_REF][START_REF] Quéléver | Effect of temperature on the formation of highly oxygenated organic molecules (HOMs) from alpha-pinene ozonolysis[END_REF][START_REF] Meusinger | Chemical and isotopic composition of secondary organic aerosol generated by α-pinene ozonolysis[END_REF][START_REF] Krechmer | Ion mobility spectrometry-mass spectrometry (IMS-MS) for on-and offline analysis of atmospheric gas and aerosol species[END_REF][START_REF] Tomaz | Structures and reactivity of peroxy radicals and dimeric products revealed by online tandem mass spectrometry[END_REF]Fang et al., 2017;Witkowski and Gierczak, 2017;Jokinen et al., 2015b;[START_REF] Nørgaard | Study of Ozone-Initiated Limonene Reaction Products by Low Temperature Plasma Ionization Mass Spectrometry[END_REF][START_REF] Bateman | Time-resolved molecular characterization of limonene/ozone aerosol using high-resolution electrospray ionization mass spectrometry[END_REF][START_REF] Walser | High-resolution mass spectrometric analysis of secondary organic aerosol produced by ozonation of limonene[END_REF][START_REF] Warscheid | Structural elucidation of monoterpene oxidation products by ion trap fragmentation using on-line atmospheric pressure chemical ionisation mass spectrometry in the negative ion mode[END_REF][START_REF] Hammes | Carboxylic acids from limonene oxidation by ozone and hydroxyl radicals: insights into mechanisms derived using a FIGAERO-CIMS[END_REF][START_REF] Kundu | High molecular weight SOA formation during limonene ozonolysis: insights from ultrahigh-resolution FT-ICR mass spectrometry characterization[END_REF]. Indeed, the use of a continuous flow reactor operating at elevated temperature, as well as a high initial concentration of reagents can induce the formation of combustion-specific products, which does not exclude their possible formation under atmospheric conditions.

Data Processing

High resolution mass spectrometry (HR-MS) generates a large amount of data that is difficult to fully analyze by sequential methods. When the study requires the processing of several thousands of molecules, the use of statistical tools and graphical representation means becomes necessary. In this study, we have chosen to use the van Krevelen diagram [START_REF] Van Krevelen | Graphical-statistical method for the study of structure and reaction processes of coal[END_REF] by adding an additional dimension, the double bond equivalent (DBE). The DBE number represents the sum of unsaturation and rings present in a chemical compound [START_REF] Melendez-Perez | A reformulated aromaticity index equation under consideration for non-aromatic and non-condensed aromatic cyclic carbonyl compounds[END_REF].

DBE = 1 + C -H/2 -O
This number is independent of the number of oxygen atoms, but changes with the number of hydrogen atoms.

Decimal values of this number, which correspond to an odd number of hydrogen atoms, were not considered in of carbon (OSc) provides a measure of the degree of oxidation of chemical compounds (alcohols, aldehydes, carboxylic acids, esters, ethers and ketones, but not peroxides) [START_REF] Kroll | Carbon oxidation state as a metric for describing the chemistry of atmospheric organic aerosol[END_REF]. This provides a framework for describing the chemistry of organic species. It is defined by the following equation:

OSc ≈ 2 O/C -H/C

Results and discussion

The oxidation of the two isomers, α-pinene and limonene (C10H16), was studied at 590 K, under atmospheric pressure, with a residence time of 2 s and a fuel concentration of 1%. Under these conditions, the formation of peroxides by autoxidation at low temperature should be optimal [START_REF] Belhadj | Experimental and kinetic modeling study of n-hexane oxidation. Detection of complex low-temperature products using high-resolution mass spectrometry[END_REF], even if the conversion of the fuels remains moderate.

Characterization of ionization sources

First, we have studied the impact of APCI and HESI sources, in positive and negative modes, on the chemical formulae detected. The HESI and APCI sources in positive and negative mode were used and their operating parameters were varied, i.e., temperature, gas flow and accelerating voltage (see Section 2). For each polarity, only ions composed of carbon, hydrogen (even numbers) and oxygen were considered. Molecular duplicates inherent to mass range overlaps were excluded. Chemical formulae with relative intensity less than 1 ppm with respect to the highest mass peak of the mass spectrum were not considered. By following these rules, we obtained a different number of ions depending on the ionization source and the polarity used. Table 1 shows the number of ions according to the experimental conditions and the discrimination rules. Each combination of ionization sources and polarity generated a set of chemical formulae. To make a meaningful comparison between the positive and negative ions data, the chemical formulae used were the precursors of the ions identified in the mass spectra. These sets have common data, but also specific chemical formulae. For a given ionization source, less than 50% of the chemical formulae are common to both polarities. In other words, between 30 and 50% of molecular species are ignored when using a single polarity. It is therefore essential to use both polarities in order to better describe all the species present. The HESI source data were compared to the APCI data (Tables 2 and3), showing an increase of the number of chemical formulae detected by 20 to 30%. This increase is characterized by a better detection of negatively ionized species and those with a higher unsaturation number (DBE). In order to evaluate further the interest for using these ionization sources, we compiled these data in Venn diagrams and proposed a visualization of these sets with a van Krevelen representation; we added the number of DBE in the third dimension. These results are presented in Tables 2 and3. In positive ionization mode, independently of the ionization source and in addition to the common molecular formulae, we detected products with an O/C ratio < 0.2 whereas in the negative ionization mode, we detected molecular formulae with an O/C ratio > 0.5. In addition to these observations, we noted that HESI is more appropriate for studying products with a large number of unsaturation (DBE > 5). Finally, for an optimal detection of the oxidation products, it is necessary to consider the transmission limits of the C-Trap. Here, we could increase by more than 60% the number of molecular formulae detected using several mass ranges for data acquisition. The most appropriate ionization polarity to be used is tight to chemical functions present in products to be detected.

We could increase by 30 to 100% the number of chemical formulae detected by using both positive and negative ionization modes. The ionization source used is also important. We could increase the number of detected chemical formulae by 20 to 30% using HESI. We believe that this approach to data validation and these results, although specific to this study, are applicable to any characterization by Orbitrap. and limonene and one which is common to both isomers. In the following, the name "only" will be used to describe the molecules specific to the oxidation of one of the terpene isomers. This common dataset represents more than 90% of the chemical formulae identified in α-pinene oxidation samples, detected with both APCI and HESI. For limonene, where the number of identified chemical formulae is greater, this common dataset ranges between 68% in APCI and 37% in HESI. In these two cases, the relatively low residence time (2 seconds) and the diversity of the chemical formulae obtained show that the oxidation of these two terpene isomers leads to the opening of the ring, a phenomenon also observed in atmospheric chemistry [START_REF] Berndt | Hydroxyl radical-induced formation of highly oxidized organic compounds[END_REF][START_REF] Kenseth | Synergistic O3 + OH oxidation pathway to extremely low-volatility dimers revealed in β-pinene secondary organic aerosol[END_REF][START_REF] Iyer | Molecular mechanism for rapid autoxidation in α-pinene ozonolysis[END_REF].

Concerning the products' molecular formulae common to both isomers, figure 1 shows that they are limited to 10 oxygen atoms. This limit is linked to α-pinene whose oxidation beyond 10 oxygen atoms remains weak (less than 2% of the totality of the identified molecules for this isomer). The autoxidation will allow keeping this similarity in the different steps of the oxidation by simple addition of two atoms of oxygen. In the case of limonene, the presence of an exocyclic double bond will increase, in a similar way to atmospheric chemistry [START_REF] Kundu | High molecular weight SOA formation during limonene ozonolysis: insights from ultrahigh-resolution FT-ICR mass spectrometry characterization[END_REF], the possibilities of oxidation and accretion. It remains however impossible, considering the size of the whole and the diversity of the isomers, to formalize all the reaction mechanisms. Nevertheless, the formation of oxidized species can be described with the help of graphical tools. The dispersion of the number of oxygen atoms per molecule shows for example that limonene oxidizes more than α-pinene (Tables 2 and3). In the case of limonene, with a HESI source, an oxygen number of up to 15 is measured, with maximum counts recorded for 10 O-atoms (Fig. 1), whereas it remains mostly less than or equal to 10 for α-pinene (Fig. 1). Moreover, this graph shows, for the products specific of limonene oxidation, a distribution centered on 11 oxygen atoms with carbon skeletons probably resulting from accretion. To verify this accretion hypothesis, we can plot the OSc as a function of the number of carbon atoms or the O/C ratio at fixed number of C-atoms (Fig. 2). One can visualize the evolution of the molecular oxidation for each carbon skeleton and the formation of oligomers. Species that are unique to one of the isomers, or common to both are differentiated using different colors. In addition to the autoxidation represented by the vertical axes for a given number of carbon atoms (Fig. 2a), we observe mechanisms of fragmentation, accretion and oligomerization between the carbon skeleton. These reaction mechanisms contribute to forming families according to the rate of oxidation and the size of their carbon skeleton. The increase in the number of oxygen atoms, but also of carbon atoms will decrease products volatility. We distinguish four families: low volatile oxygenated organic aerosols (LV-OOA), semi-volatile oxygenated organic aerosols (SV-OOA), biomass burning organic aerosols (BBOA) and water-soluble organic carbons (WSOC) following a classification proposed in the literature [START_REF] Kroll | Carbon oxidation state as a metric for describing the chemistry of atmospheric organic aerosol[END_REF].

If
we analyze the two sets of molecules from the APCI and HESI sources, positive and negative ionization modes combined, we find that nearly 73% of the molecules are linked to each other by a single difference of 2 oxygen atoms which reflects an autoxidation mechanism (Fig. 3). We can measure the amplitude of autoxidation for each carbon skeleton from the OSc vs. O/C space. For the two isomers, for which the initial carbon number (C) is equal to 10, one can observe (Fig. 2b) two autoxidation routes with an even and odd number of oxygen atoms. This parity distinction is initially present for the two main radicals, ROO˙ and RO˙, responsible for autoxidation. However, the termination and propagation reactions will change the oxygen parity to form a new distribution where the parity links between radicals and molecules are lost, making any interpretation of radical oxidation route impossible (Fig. 3). The autoxidation mechanisms indicated by arrows in Figure 2b characterize systematic hydroperoxy chain terminations which do not change oxygen parity. HESI data showed an equivalent distribution of the two oxygen parities in molecular products (odd: 51%; even 49%) therefore confirming a lack of selectivity of the reaction mechanisms with respect to the parity of the radicals. Chain propagation

RO • +RH ROH + R • ROO • + RH ROOH + R • ROO • QOOH • KHP + OH • Chain termination R • + RO • ROR ROO • + R'OO • ROOR' / R

Combustion versus atmospheric oxidation

We have explored potential chemical pathways related to autoxidation in the previous section. For this purpose, we have performed experiments under cool flame conditions (590 K). This autoxidation mechanism is also present in atmospheric chemistry, but it is only recently that it has been found that this mechanism could be one of the main formation pathways for SOA [START_REF] Savee | Carbon radicals. Direct observation and kinetics of a hydroperoxyalkyl radical (QOOH)[END_REF]. Studies have described this mechanism in the case of atmospheric chemistry with the identification of radicals and molecular species [START_REF] Tomaz | Structures and reactivity of peroxy radicals and dimeric products revealed by online tandem mass spectrometry[END_REF]. However, previous studies of the propagation of this reaction mechanism are mainly focused on the initial skeletons of the 10-carbon isomers, whereas the other carbon skeletons are also concerned by autoxidation. It is therefore useful to evaluate the proportion of products of autoxidation among the total species formed.

We propose a new approach which consists in assessing a set of molecules mainly resulting from autoxidation against different sets of experimental studies related to atmospheric chemistry. The objective is to compare the similarity between these conditions on the basis of autoxidation by considering different experimental parameters chosen for their diversity. For this purpose, we selected, a HESI ionization source, better adapted to the electronegativity of the oxidized molecules, as well as to higher m/z. Moreover, we have already shown in this paper that the HESI (+/-) ionization source is better suited for detecting autoxidation products (detection of 96% of the total chemical formulae observed in autoxidation).

Among published atmospheric chemistry studies of terpenes oxidation, we have selected 15 studies, 4 for α-pinene and 11 for limonene. The data were acquired using different experimental procedures (methods of oxidation, techniques of characterization). Table 4 summarizes all the experimental parameters related to the selected studies. The data are from the articles or files provided in the supporting information S2. In these studies, oxidation was performed only by ozonolysis with different experimental conditions that gather the main methods described in the literature: ozonolysis, dark ozonolysis, ozonolysis with OH scavenger, ozonolysis with or without seed particles. We considered that the methods of analysis by mass spectrometry did not modify the nature of the chemical species but only their relative importance, because of the type of ionization and the sensitivity of the instruments. The combination of data obtained using (+/-) HESI gives a rather complete picture of the autoxidation products.

First, we compared the data from ozonolysis studies of each terpene and identified similarities through the Venn diagram. For studies with two ionization sources, duplicate chemical formulae were removed. We selected the three most representative studies, by the number of the chemical formulae detected. Then, we compared the set of chemical formulae identified after ozonolysis to those produced in low-temperature combustion, the objective being (i) to highlight similarities in terms of products generated by the two oxidation modes and (ii) to identify data resulting from autoxidation.

For α-pinene oxidation, the four studies identified 566 chemical formulae, all polarities combined. Only the study by [START_REF] Meusinger | Chemical and isotopic composition of secondary organic aerosol generated by α-pinene ozonolysis[END_REF]was performed in positive mode and none of the studies reported data obtained with two ionization modes (+/-). For limonene oxidation, the three studies identified 1862 chemical formulae. Only the studies by [START_REF] Walser | High-resolution mass spectrometric analysis of secondary organic aerosol produced by ozonation of limonene[END_REF][START_REF] Bateman | Time-resolved molecular characterization of limonene/ozone aerosol using high-resolution electrospray ionization mass spectrometry[END_REF] used (+) and (-) ionization modes. In the case of limonene, oxidation for which accretion is more important than for α-pinene, and for which a greater number of chemical formulae were identified, the similarities is more important. These results are presented in Figure 4 where the ionization polarity used in each study is specified. (see conditions in Table 1). Each digit of the binary numbers in parentheses identifies the datasets being compared.

For a-pinene, no global similarity is observed for the detected chemical formulae. Different hypotheses can be offered to explain this result. Among them, the number of chemical formulae identified per study remains limited (a few dozen to several hundred) and these small datasets are sometimes restricted to specific mass ranges (e.g. For limonene, 86 chemical formulae are common to the three studies considered here. In this dataset, nearly 43% of the chemical formulae present a relation similar to that of autoxidation (simple difference of two oxygen atoms)

with less than 5% of fragments (Nb C<10). This result seems to indicate that autoxidation dominates.

One can then ask if reaction mechanisms common to atmospheric and combustion chemistry can generate, despite radically different experimental conditions, a set of common chemical formulae and if in this common dataset, a common link, characteristic of autoxidation, is observable?

To address that question, we compared all the previous results, for each of these terpenes to those obtained under the present combustion study. The comparisons were made using the HESI. One should remember that the oxidation conditions in a JSR were chosen in order to maximize low-temperature autoxidation. Again, we used Venn diagrams to analyze these datasets composed of 1590 chemical formulae in the case of α-pinene and 2857 chemical formulae in the case of limonene. The results of these analyses are presented in Figure 5. 1). specify again that different reaction mechanisms can cause the observed similarities. However, the preponderance of autoxidation in so-called cold flame combustion is obvious, and in atmospheric chemistry, this reaction mechanism remains competitive or dominates [START_REF] Crounse | Autoxidation of organic compounds in the atmosphere[END_REF]Jokinen et al., 2014b). If we search for an autoxidation link between these 86 chemical formulae, we observe that 43% of these chemical formulae meet this condition: difference of two oxygen atoms, at constant number of carbon and hydrogen atoms. More precisely, these molecules are centered in a van Krevelen diagram on the ratios O/C=0.6 and H/C=1.6, in the range 0.25 < O/C < 0.8 and 1.25 < H/C < 1.85. All oxidized molecules associated with this dataset are presented in Figure 6.

The dispersion of the chemical formulae, far from being random, remains consistent with an autoxidation mechanism where the number of carbon atoms is constant. Specificities and similarities of these two oxidation modes were further investigated by plotting the distribution of the number of oxygen atoms in detected chemical formulae (Fig. 7). Indeed, the distribution of the number of oxygen atoms allows, in addition to the Van Krevelen diagram, to provide some additional details on these two modes of oxidation. It is in ozonolysis that we observe the chemical formulae having the largest number of oxygen atoms. There, oxidation proceeds over a long reaction time where the phenomenon of aging appears by promoting accretion or oligomerization. One should keep in mind that a short residence time (2 s) was used in JSR experiments. A short residence time should limit accretion which is favored under simulated atmospheric oxidation experiments performed over long reaction times. In combustion, the number of oxygen atoms remains limited to 18, with a lower number of detected chemical formulae compared to the case of ozonolysis. However, it is in combustion that we observe the highest O/C ratios, indicating the formation of the most oxidized products. This difference, however, does not affect the similarities between the chemical formulae detected in the two modes of oxidation. Finally, the analysis of the parities in oxygen atoms, very similar to the three datasets, confirms that the reaction mechanisms presented in Figure 3 do not allow a simple link to be established between the parity of radicals and that of molecular products. chemical relationship between molecules). We selected existing databases on ozonolysis with a purposive diversity of experimental conditions and compared these data to those obtained in low-temperature combustion conditions where autoxidation is important (cool flame combustion at 590 K). We showed that a significant portion of the chemical formulae were common to both atmospheric and combustion products. Surprisingly, the numerous oxidation mechanisms and the isomerization of chemical species proceeding under these two conditions did not lead to diverging data, but, on the contrary, to similarities. More than 35% of the chemical formulae detected in combustion chemistry experiments using a JSR have been detected in the studies carried out under atmospheric conditions. Finally, we have outlined the existence of a substantial common dataset of autoxidation products. This result tends to show that autoxidation is indeed inducing similarity between atmospheric and combustion products.

  https://doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.

  this study. Then, duplications of chemical formulae in the O/C vs. H/C space are eliminated. The oxidation state https://doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.

  https://doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.

  https://doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.

  doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License. Table 3. Representation of the mass spectrometry data characterizing the oxidation of α-pinene and 239 limonene (ionization source: HESI positive and negative modes, JSR experiments). 240 products detected in a JSR 242 In order to compare the oxidation of α-pinene and limonene, we compiled the positive and negative ionization data 243 obtained with APCI (Table2) and HESI (Table3) ionization sources to obtain a more exhaustive database. For 244 the APCI and HESI sources, we distinguished three datasets, two of which are specific to the oxidation of α-pinene /doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.

Figure 1 :

 1 Figure 1: Distribution of α-pinene and limonene autoxidation products as a function of their oxygen content (ionization source: HESI, combined positive and negative modes data).

  doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.

Figure 2 :

 2 Figure 2: Overview of the distribution of limonene and α-pinene oxidation products observed in a JSR: (a) OSc versus carbon number in detected chemical formulae from APCI and HESI data. (b) OSc versus O/C atomic ratio for a carbon number of 10; index of the products: number of hydrogen atoms. Arrows indicate autoxidation from a C6H16 isomer, according to the oxygen parity in products.

  https://doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License. of the oxygen atom number in a reaction mechanism Even (49%) and odd (51%) number of oxygen atoms RH ROO •

Figure 3 :

 3 Figure 3: Autoxidation reaction mechanisms in combustion (left) and in the atmosphere (left and right). * indicates a change in number of oxygen atoms.

  https://doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License. Table 4. Experimental settings of 15 oxidation studies of two terpenes under atmospheric conditions and cool doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.

Figure 4 :

 4 Figure 4: Venn diagrams for comparing the oxidation results from ozonolysis of (a) α-pinene and (b) limonene

  doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.[START_REF] Quéléver | Effect of temperature on the formation of highly oxygenated organic molecules (HOMs) from alpha-pinene ozonolysis[END_REF]: C10 to C20). In the case of studies carried out with an NO3 -source, sensitive to HOMS, produced preferentially by autoxidation, we note that nearly 50% of the chemical formulae (10/22) are linked by a simple difference of 2 oxygen atoms.

  It turned out that 301 chemical formulae were identified for α-pinene and 735 chemical formulae for limonene to be common to oxidation by ozonolysis and combustion. This represents 31% of the chemical formulae for the ozonolysis of α-pinene and 31% for those of limonene ozonolysis. For α-pinene, the similarities compared to combustion are specific to each study:[START_REF] Deng | Temperature and acidity dependence of secondary organic aerosol formation from α-pinene ozonolysis with a compact chamber system[END_REF] 69%[START_REF] Meusinger | Chemical and isotopic composition of secondary organic aerosol generated by α-pinene ozonolysis[END_REF] 46%[START_REF] Quéléver | Effect of temperature on the formation of highly oxygenated organic molecules (HOMs) from alpha-pinene ozonolysis[END_REF]7%[START_REF] Krechmer | Ion mobility spectrometry-mass spectrometry (IMS-MS) for on-and offline analysis of atmospheric gas and aerosol species[END_REF] 23%. Chemical formulae common to all studies were not identified. This lack of similarity may be due to a partial characterization of the chemical formulae, a weaker oxidation of α-pinene with an ionization mode less favorable to low molecular weights.

Figure 5 :

 5 Figure 5: Venn diagrams comparing the oxidation results from ozonolysis and combustion of (a) α-pinene and (b) limonene (see conditions in Table1).

  103

  doi.org/10.5194/acp-2021-964 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License. For limonene, the similarities with combustion are more important and less spread out. They represent for the different studies: (Kundu et al., 2012) 65% (Walser et al., 2008) 88% (Tomaz et al., 2021) 81%. Moreover, there is a common dataset of 86 chemical formulae which can derive from autoxidation mechanisms. It is necessary to

Figure 6 :

 6 Figure 6: Van Krevelen diagram showing specific and common chemical formulae detected after to oxidation of limonene by ozonolysis and combustion.A 3-D representation of all limonene oxidation data is given in Fig.S1(Supplement) where DBE is used as third dimension. From that figure, one can note that products with higher DBE are preferably formed under JSR conditions, i.e. at elevated temperature. The corresponding chemical formulae could correspond to carbonyls and / or cyclic ethers ( • QOOH → carbonyl + alkene + OH and / or cyclic ether + OH).

Figure 7 :

 7 Figure 7:Oxygen number distribution for all the molecules identified for the oxidation of limonene: only in combustion, only in ozonolysis and common to both processes. For each of these distributions, the parity in oxygen number is specified5 ConclusionWe have analyzed and compared, thanks to different mathematical and visualization tools, the oxidation of αpinene and limonene under simulated tropospheric and low-temperature oxidation conditions. We have restricted the field of this study to organic compounds (CnHyOz) in order to study specifically autoxidation. This work is in the continuity of recently published studies which established the importance of autoxidation under tropospheric oxidation and low-temperature combustion conditions. Faced with the analytical limits of the chemical speciation of thousands of molecules, and in order to complete the work carried out on the reaction mechanisms of the first steps of oxidation, we have developed a global approach based on the study of families of chemical compounds present in these large sets (van Krevelen diagram, oxygen number distribution, oxidation state of carbon -OSc,

Table 1 .

 1 Number of ions detected for each source in positive and negative modes.

	Ionization source	α-Pinene		Limonene
	APCI	646 (+)	503(-)	1321(+)	1346(-)
	HESI	594(+)	693(-)	1017(+)	1864(-)

Table 2 .

 2 Representation of the mass spectrometry data characterizing the oxidation of α-pinene and 235 limonene (ionization source: APCI positive and negative mode, JSR experiments).

	236
				α-Pinene			Limonene
	Source and mode	APCI	+		APCI -	APCI +		APCI -
	Number of	646			503	1321		1346
	compounds
	Distribution	specific to		common	specific to	specific to	common	specific to
		positive mode	319	negative mode	positive mode	805	negative
		327			174	516		mode
								542
	Venn graph

Acknowledgements

The authors gratefully acknowledge funding from the Labex Caprysses (ANR-11-LABX-0006-01), the Labex Voltaire (ANR-10-LABX-100-01), CPER, and EFRD (PROMESTOCK and APPROPOR-e projects). We also thank Matthieu RIVA for sharing his experimental data on limonene oxidation.


B Aubouin-Pairault
email: bob.aubouin-pairault@gipsa-lab.fr

A Perodou
email: arthur.perodou@ec-lyon.fr

C Combastel

A Zolghadri
email: ali.zolghadri@ims-bordeaux.fr.

Resilient tube-based MPC for Cyber-Physical Systems Under DoS Attacks 1

Keywords: Cyber-Physical Systems (CPS), Denial of Service (DoS) attack, Resilience, Robustness, Model Predictive Control (MPC), Set-membership, Zonotopes

This paper proposes a resilient and robust model predictive control (MPC) scheme for a class of Cyber-Physical Systems (CPS) subject to state and input constraints, unknown but bounded disturbances and Denial of Service (DoS) attacks. The attacker blocks the controller to actuator communication and the attacks are assumed to be time-limited. The control is designed by extending a robust tube-based MPC, where a new type of invariant set, namely µ-step Robust Positively Invariant (µ-RPI) set, is introduced to deal with resilience. A set-based method is then developed for the control scheme to ensure both resilience to DoS attacks while preserving robustness to bounded disturbances. A computational algorithm is derived and a numerical example is provided to illustrate the potential of the proposed approach.

INTRODUCTION

Cyber-Physical Systems (CPSs) are next-generation engineered systems with deep integration of computation, communication and networking, physical processes, and control systems [START_REF] Lee | The past, present and future of cyberphysical systems: A focus on models[END_REF], [START_REF] Poovendran | Special issue on cyber-physical systems [scanning the issue[END_REF], [START_REF] Allgöwer | Position paper on the challenges posed by modern applications to cyber-physical systems theory[END_REF]). Varying widely in complexity and scale, CPSs concern many technological areas, including aerospace, automotive, energy, chemical industry, transportation, or healthcare. While safety and security have been traditionally addressed separately, some integrated approaches have been recently reported (see for instance [START_REF] Ji | Harmonizing safety and security risk analysis and prevention in cyber-physical systems[END_REF] and the references therein). In this perspective, this paper focuses on mitigating the impact of denial of services (DoS) attacks on systems evolving in an uncertain environment. In the literature, one can find many investigations for mitigation of DoS attacks. Among others, in [START_REF] Gupta | Dynamic games with asymmetric information and resource constrained players with applications to security of cyberphysical systems[END_REF] a method is proposed based on game theory where both attacker and controller are modeled as players of a game. De [START_REF] De Persis | Input-to-state stabilizing control under denial-of-service[END_REF] used event-triggered control, and the approach reported in [START_REF] Amin | Safe and secure networked control systems under denial-ofservice attacks[END_REF] takes advantage of optimal control. Nevertheless few of them deal with systems subject to state and input constraints. [START_REF] Sun | Resilient model predictive control of cyber-physical systems under DoS attacks[END_REF] proposed a resilient MPC to address DoS attacks for a constrained system but without taking disturbances into account. In this paper, a similar resilient MPC is extended to disturbed systems with the aim of ensuring a joint resilience and robustness goal. By resilience, it is meant the system's ability to contain the maximal impact of anomalies, such as attacks, and to recover to an acceptable performance level.

In order to characterize the impact of uncertainties, a set-membership approach is used. This deterministic approach has been used in various ways when dealing with systems subject to unknown but bounded uncertainties. For instance, state bounding observation through Zonotopic Kalman Filters (ZKF) is considered in [START_REF] Combastel | Zonotopes and Kalman observers: Gain optimality under distinct uncertainty paradigms and robust convergence[END_REF]. In [START_REF] Mayne | Robust output feedback model predictive control of constrained linear systems[END_REF] a robust model predictive controller is proposed, and in [START_REF] Le | Robust tube-based constrained predictive control via zonotopic set-membership estimation[END_REF] a zonotopic tube-based approach is used to control a system subject to disturbances and measurement noise. In the last two papers, a particular class of invariant sets, the Robust Positively Invariant (RPI) set, is defined to ensure the robust convergence of the control scheme in the presence of bounded uncertainties. [START_REF] Franze | Resilient model predictive control for constrained cyber-physical systems subject to severe attacks on the communication channels[END_REF] also studied set membership theory and MPC to control constrained and disturbed linear systems subject to different possible attacks. The resilient approach we propose mainly relies on zonotopic rather than ellipsoidal sets. It ensures that the state trajectory comes back within some set after a DoS attack sequence, as in [START_REF] Franze | Resilient model predictive control for constrained cyber-physical systems subject to severe attacks on the communication channels[END_REF]. Moreover, an explicit bound for the state trajectory during attacks is also considered.

In this paper, a new kind of sets, namely µ-step Robust Positively Invariant (µ-RPI) sets, is introduced to address the resilience issue, that is, ensuring that an exit from nominal operation induced by the possibly repeated occurrence of attacks remains limited in time and/or magnitude. Then, a resilient tube-based model predictive control for uncertain linear discrete-time systems subject to DoS attack and bounded disturbances is proposed. The resulting MPC scheme relies on the off-line computation of such an µ-RPI set for which an algorithmic solution is provided, taking into account all design constraints. The paper is organized as follows. After some preliminaries given in Section 2, the problem is stated in Section 3. Section 4 describes the proposed control scheme and Section 5 presents a method to compute the required µ-RPI set. Finally an illustrative example is provided in Section 6, and Section 7 provides some concluding remarks.

PRELIMINARIES

Notations:

The symbols R, N and N + represent respectively the sets of real numbers, natural numbers, and positive integers. k 1 , k 2 denotes the set of integers between, and including, k 1 and k 2 . The operators ⊕ and ⊖ denote the Minkowski sum and difference. Given the sequence of sets {S i , ⊂ R n } b i=a , the notation b i=a S i = S a ⊕ ... ⊕ S b is used. P ≻ 0 (resp. P ⪰ 0) denote a positive definite (resp. semi-definite) matrix. ρ(A) denotes the spectral radius of the matrix A i.e. the largest eigenvalue in absolute value. Given a vector x ∈ R n and P ≻ 0, x ⊤ denotes the transpose of the vector and the P -weighted norm is ||x|| P =

√

x ⊤ P x. The notations x k1|k2 and u k1|k2 are used to describe state and control input predicted for time k 1 and calculated at time k 2 .

A zonotope ⟨c, R⟩ ⊂ R n with center c ∈ R n and generator matrix R ∈ R n×p is a polytopic set defined as the affine image of the unit hypercube

[-1, +1] p ⊂ R p by R: ⟨c, R⟩ = {c + Rs, ||s|| ∞ ≤ 1}. The Minkowski sum of two zonotopes ⟨c 1 , R 1 ⟩ and ⟨c 2 , R 2 ⟩ is the zonotope ⟨c 1 , R 1 ⟩ ⊕ ⟨c 2 , R 2 ⟩ = ⟨c 1 + c 2 , [R 1 , R 2 ]⟩.
The linear image of ⟨c, R⟩ by L is the zonotope L ⊙ ⟨c, R⟩ = ⟨Lc, LR⟩. The interval/box hull of ⟨c, R⟩ is the interval c ± |R|1 where c ± r denotes the interval [c -r, c + r], |.| is the elementby-element absolute value operator, and 1 is a column of ones with appropriate size, [START_REF] Combastel | A state bounding observer based on zonotopes[END_REF].

PROBLEM STATEMENT

System Dynamics

A linear dynamic subject to bounded disturbances is first considered:

x k+1 = Ax k + Bu k + Dw k , (1) where x k ∈ R n , u k ∈ R m ,
and w k ∈ R l are respectively the system state, the control input and the disturbances evaluated at time k. It is assumed that the disturbances are bounded: w k ∈ W and the system state and control input satisfy constraints x k ∈ X and u k ∈ U, where X and U are compact sets containing the origin as an interior point. Assumption 1. There exists a state-feedback gain K such that A + BK is stable, that is ρ(A + BK) < 1.

Attack Model

The equation ( 1) is modelling the attack-free dynamic of a CPS where actuator and controller are spatially separated (e.g. remote control) while being connected through communication channel. In this paper, the case of a DoS attack occurring on this channel is considered. More precisely, the scenario considered here is time limited attacks which block the controller-actuator (C-A) communication channel. The indicator variable ν k , k ∈ N, is defined as follows:

ν k =
1 while there is no attack, 0 when an attack occurs.

(2)

Then, the applied input is u k = ν k u c k where u c k is the control input computed by the controller. Assumption 2. ∀k j ∈ N, the constraint on the occurrence of attacks is the following:

kj +N -1 k=kj (1 -ν k ) ≤ M < N, (3)
where M denotes the maximal number of sample attack instances on the time interval from k j to k j + N -1.

This assumption could model power-limited attacks where the attacker can launch a maximal occurrence of attacks over a time period. M consecutive attacks are not necessarily the worst-case scenario, the assumption considers all the possible scenarios with at most M attacks over N sample times. Indeed, under assumption 2, distinct consecutive DoS attack sequences are possible within the considered time horizon, as illustrated in Fig. 3. In this paper, the number N matches with the time-horizon of the MPC scheme that will be introduced in next section.

Resilient sets and problem statement

In order to deal with the resilience property for a system subject to disturbances and attacks, a new kind of invariant set is proposed in this subsection. First, the definition of a Robust Positively Invariant (RPI) set is recalled: Definition 1. (RPI set, [START_REF] Rakovic | Invariant approximations of the minimal robust positively invariant set[END_REF]). Given a dynamic system x k+1 = f (x k , w k ), the set S is said to be Robust Positively Invariant (RPI), if f (x, w) ∈ S for all x ∈ S and all w ∈ W.

A resilient and robust invariant set, that both generalizes the usual RPI set and the µ-step invariant set originally proposed in [START_REF] Sun | Resilient model predictive control of cyber-physical systems under DoS attacks[END_REF], is now introduced: Definition 2. (µ-RPI set). Given a system Σ modelled as

x k+1 = f (x k , w k ), if ∃µ ∈ N + such that the implication (4) is satisfied, then the set S is said to be µ-step Robust Positively Invariant (µ-RPI set) for Σ, ((x 0 ∈ S) ∧ (∀k ∈ N, w k ∈ W)) ⇒ ∀k ∈ N, x µ+k ∈ S. (4)
Given the system (1) subject to a DoS attack modeled by (2), the problem addressed is that of finding a resilient and robust model predictive control scheme for which an µ-RPI set exists, and compute this set. By this way, the purpose is to obtain a provenly stable control scheme which is jointly resilient to attacks as described in the paragraph 3.2 and robust to any arbitrary bounded disturbances w k ∈ W.

A ROBUST AND RESILIENT MPC SCHEME

In this section, a standard robust tube-based MPC is first briefly recalled. Then, a recent proposition of MPC scheme providing resilience to a non-disturbed LTI system under DoS attacks is reviewed. Finally, an original MPC scheme combining both robustness and resilience is proposed.

Tube-based MPC

In standard robust tube-based MPC, the objective is to control the system (1) such that its states remain in a tube centered on a nominal trajectory. Consider the following nominal model:

xk+1 = Ax k + B ūk . (5)
Then, the general idea is to ensure that the disturbed system state x k stays close to xk , for any arbitrary disturbances belonging to W. To achieve this goal, the nominal system is controlled using MPC with tightened constraints and the control law ( 6) is applied to the disturbed system:

u k = ūk + K(x k -xk ). ( 6
)
If K is chosen such that (A + BK) is stable, the difference x k -xk can be bounded in a Robust Positively Invariant (RPI) set and thus creates a tube around the nominal state trajectory that encloses the disturbed state trajectories as described in [START_REF] Rawlings | Model Predictive Control Theory and Design[END_REF]. Note that this is a simple version of the tube-based MPC where the nominal system is completely independent of the state trajectory and that improved versions of this control method exist (see [START_REF] Rawlings | Model Predictive Control Theory and Design[END_REF]).

Resilient MPC

This paragraph is intended to provide a short synopsis of the resilient MPC for non-disturbed systems as proposed in [START_REF] Sun | Resilient model predictive control of cyber-physical systems under DoS attacks[END_REF]. In a classical MPC scheme, at each step time, the optimal sequence of control inputs u * k is computed over the whole N samples prediction horizon and only the first term u * k|k is applied to the system. In order to mitigate the impact of DoS attacks, [START_REF] Sun | Resilient model predictive control of cyber-physical systems under DoS attacks[END_REF] made the assumption that the actuator side is able to store data and the computed control sequence {u * k|k , u * k+1|k , .., u * k+N -1|k } is sent to the actuator at each sample time. Thus, when a DoS attack impedes the reception of the optimal control at some time k, the actuator can still apply the next control stored in the previously received control sequence. Denoting k j the last sample time before an attack occurs, the control applied to the disturbed system is u * k|kj which leads to the dynamic:

x k+1 = Ax k + Bu * k|kj . (7)
In [START_REF] Sun | Resilient model predictive control of cyber-physical systems under DoS attacks[END_REF], the stability of ( 7) is proved under assumption 2 and for a region of attraction X M N . The later differs from that of a classical MPC by taking into account the case of an attack happening at a time when no control sequence is stored in the actuator side (just after the initialization of the controller, for instance). This set is such that for all x 0 ∈ X M N , there exists a control sequence u M = {0, .., 0

M , u M , .., u N -1 N -M } satisfying x N ∈ X f under the constraints u i ∈ U, ∀i ∈ M, N -1 and x i ∈ X, ∀i ∈ 1, N -1 ,
where X f denotes the terminal set of the MPC problem.

Toward a robust and resilient control law

To provide both resilience and robustness, a new control law combining the previous tube-based MPC ( §4.1) and resilient MPC ( §4.2) is proposed. First, the constrained optimization problem which is solved to control the nominal system (5) is introduced:

MPC optimization: The cost function is given by:

V N (x k , ūk ) = ||x k+N |k || 2 S + N -1 i=0 ||x k+i|k || 2 Q + ||ū k+i|k || 2 R
where ūk = {ū k|k , ūk+1|k , .., ūk+N-1|k }. and Q, R, S ≻ 0.

The constraints applied to the state and the control input of the nominal system are X = X ⊖ Z + and Ū = U ⊖ KZ + , where Z + denotes a set bounding the difference z k = x kxk . Those tightened constraints applied to the nominal system ensure that the disturbed system (1) meets the specification x k ∈ X and u k ∈ U.

In order to ensure the stability and feasibility of the considered MPC problem, a terminal set Xf is also introduced (see [START_REF] Rawlings | Model Predictive Control Theory and Design[END_REF] for details). The terminal cost ||..|| S and the terminal set Xf are chosen consistently according to assumption 3: Assumption 3. There exists a state-feedback gain K such that, for the system xk+1 = (A + BK)x k :

• There exists a terminal set Xf ⊂ X also satisfying

(A + BK) Xf ⊂ Xf and K Xf ⊂ Ū, • ||(A + BK)x k || 2 S + ||x k || 2 Q + ||K xk || 2 R ≤ ||x k || 2 S .
Then, the MPC constrained optimization is given by: ū

* k = min ūk V N (x k , ūk ) s.t. xk|k = xk , xk+i+1|k = Ax k+i|k + B ūk+i|k , ūk+i|k ∈ Ū for i ∈ 0 : N -1 , xk+i|k ∈ X for i ∈ 1 : N -1 , xk+N|k ∈ Xf . (8)
Resilience: Since attacks can occur, the robust control (6) may not be received by the actuators at some time steps. Similarly to [START_REF] Sun | Resilient model predictive control of cyber-physical systems under DoS attacks[END_REF] and [START_REF] Franze | Resilient model predictive control for constrained cyber-physical systems subject to severe attacks on the communication channels[END_REF], it is assumed that the actuator side is able to store data. Then, since only past measured states are available at each sample time k, the sequence {ū * k|k + K(x k -xk ), ū * k+1|k , .., ū * k+N -1|k } is sent to the actuator at each sample time k. Denoting k j the last sample time before an attack occurs and using ν k as defined in (2), the control applied to the disturbed system is:

u k = ū * k|kj + ν k K(x k -xk )
For the sake of simplicity, x k is assumed available. Otherwise, an observer should be used. Then, the robustness to the related observation error should be properly handled. Moreover, in a similar fashion to [START_REF] Amin | Safe and secure networked control systems under denial-ofservice attacks[END_REF]), the communication protocol is assumed to be acknowledgment-based like e.g. the TCP protocol [START_REF] Kumar | Survey on transport layer protocols: TCP & UDP[END_REF]). Then, the controller has access to k j and the nominal system can be controlled using ūk = ū * k|kj . With this control, the nominal system follows the dynamic (7) already studied by [START_REF] Sun | Resilient model predictive control of cyber-physical systems under DoS attacks[END_REF]. The controlled nominal system is thus stable for a region of attraction XM N . This region of attraction XM N is the same as in subsection 4.2, except that that the constraints on the states and the control inputs are given by the tightened constraints Xf , Ū and X.

Robustness and Resilience: The proposed approach relies on a dedicated management of the state error between the disturbed (1) and nominal (5) system models, z k = x k -xk , which follows the dynamics:

z k+1 = (A + ν k BK)z k + Dw k (
9) This is the subject of the main results developed in the next section.

MAIN RESULTS

In this section, an algorithmic solution is provided to test if a given RPI set for the dynamic (10),

z k+1 = (A + BK)z k + Dw k , (10
) is also an N -RPI set for the dynamic (9) under attack scenarios characterized by the pair (M, N ). As a byproduct, it will thus be possible to determine the maximal value of M < N ensuring a resilient MPC control. Moreover, a solution for computing a set Z + bounding the trajectories of the state error z k is proposed. Z + is represented as an intersection of halfspaces and characterizes the worstcase impact of the specified disturbances and attacks with respect to the nominal state trajectory.

Testing if an RPI set is also an N -RPI set

In this section, P denotes an RPI set for the system (10). It can be computed as an approximation of the minimal RPI set of (10), as proposed in [START_REF] Rakovic | Invariant approximations of the minimal robust positively invariant set[END_REF].

The aim of this paragraph 5.1 is to test if P is also an N-RPI set for system (9). Indeed, since the condition (3) applies to sequences of N steps, a solution to ensure the resilience to repeating sequences of any length consists in focusing on µ-RPI sets with µ ≤ N (thus, on N -RPI sets), as illustrated in Fig. 1.

Fig. 1. Illustration of the usage of the N -RPI

To achieve this, the proposed method consists in testing if, for each considered attack scenario, the trajectory gets back into P in a maximum of N steps. Instead of testing all the possible attack scenarios, a sufficient subset is introduced in the sequel. In addition to reduce the algorithmic complexity, there is no loss of generality to prove, or disprove, that P is an N -RPI set when testing only this subset. This relevant subset is obtained as follows: Using the above rules, a sufficient subset of scenarios to be tested can be constructed. First, starting from an RPI set P at k = 0 and with one sample attack happening at this time, a minimal number k b of steps (up to N ) required for the system trajectories to provenly get back to P is determined. Then, the expression of the reachable set is computed as follows:

   Z 1 0 = P, Z 1 k+1 = AZ 1 k ⊕ DW, if ν k = 0, Z 1 k+1 = (A + BK)Z 1 k ⊕ DW, otherwise. ( 11
)
Where Z j k is the set reached during the k th step for a scenario j. To test the scenarios with two sample attack instances, the initial set is given by Z 2 0 = Z 1 1 and the scenarios considered are the ones where an attack occurs at time k with k ∈ 0, k b -2 . The number of sample attacks is then incremented until it reaches M which is necessary for the test to succeed unless a set-inclusion test fails.

Algorithm 1 is a recursive function that tests if an RPI set P is also an N -RPI set for attack scenarios described by the integers (M, N ). Each call of the recursive function tests all the required scenarios, starting from the set Z 0 with an attack happening between step 0 and step n step -1, where n step is a number of steps passed as input argument (with initial value one) when calling the isNRPI function. For a scenario with an attack happening at time k a , and the trajectories entering in the RPI set in less than N steps, in k b steps for instance. The recursive function is recalled to test new scenarios starting from the set S = Z ka+1 and it tests that the trajectories will come back in the RPI set in less than N -k a -1 steps with an attack happening between the steps 0 and k b -2-k a . The function stops if a trajectory does not enter in the RPI before the maximum allowed step number or if there is no more attack to test i.e. the last scenarios tested already include the cases with M sample attacks.

For the sake of illustration, the operation of Algorithm 1 is exemplified in Fig. 2 where:

• During the first call, a scenario starting from P with an attack at the first step is considered. The trajectories get back to P in 2 steps. • The function is recalled and a scenario starting from Z 1 1 with an attack at the first step is considered. The trajectories get back to P in 3 steps.

• During the third call of the function, two scenarios are studied. Both start from the set Z 2 1 and one contains an attack at the first step and the other one an attack at the second step. The two trajectories from the scenarios come back in P in less than N -2 steps which is the limit for those trajectories.

Algorithm 1 Test if a RPI P is a N -RPI for a given number M of sample attacks over the horizon N

1: function [test, b] = isNRPI(M , N , P, Z 0 , n step , b) 2: if (Z 0 , n step , b) is undefined then ▷ Initialization 3:
Z 0 ← P ▷ where P has to be a RPI set 4:

n step ← 1 ▷ First attack happens at time k=0 5: b = (-∞)1 ▷ col. vector with -∞ elements 6: end if 7: if (M = 0) then ▷ Stop criterion for recursive calls 8: test ← T rue 9: else 10: k a ← 0 ▷ Sample attack time-step 11:
Cond 1 ← T rue 12:

while Cond 1 do 13:

i ← 0 ▷ Time-step index end if

23: b ← U pdateBound(H, b, Z) ▷ Update Z+ 24: i ← i + 1 25: Cond 2 ← (i < N ) ∧ ¬(Z ⊂ P) 26:
end while 27:

k b ← i 28:
if (Z ⊂ P) then 29:

▷ When back to P in less than N steps, then 30:

▷ recursive call with one less sample attack 31:

▷ and k a less steps to get back to P starting 32:

▷ from S:

33: (test, b) ← isNRPI(M -1, N -k a -1, P, S, k b -1 -k a ,
k a ← k a + 1 38: Cond 1 ← (k a < n step ) ∧ (test) 39:
end while 40: end if 41: return test, b

The implementation of Algorithm 1 makes use of zonotopes to represent the RPI set and the sets reachable by the trajectories. This representation is well-suited in this case since it is closed for the Minkowski sum and affine image operators which can be easily computed (see §2). To test the containment of zonotopes at step 25 of Algorithm 1, the algorithm presented in [START_REF] Kulmburg | On the co-NPcompleteness of the zonotope containment problem[END_REF] and implemented in the CORA 2021 toolbox is used. Note also that the step 23 in Algorithm 1 is used to compute a set bounding all the trajectories as further described in the next paragraph 5.2.

Bound of trajectories

The aim of this paragraph is to compute an outer approximation Z+ of the set Z + bounding the trajectory error with respect to the nominal model, under both bounded disturbances and the specified attacks. Following (9) with z 0 ∈ P, where P is N -RPI as tested by Algorithm 1, Z + is given by ( 12):

Z + = n Z i=1 Z i , (12)
where Z i , i ∈ 1, n Z denotes the collection of the sets Z computed at steps 18 and 21 during all the iterations within Algorithm 1 i.e. during the recursive exploration of the required scenarios. Since zonotopes are not closed for union, Z + is not a zonotope. To compute an outer approximation Z+ of Z + , a polytope in H-representation (i.e. intersection of half-spaces) can be used, Z+ = {z | Hz ≤ b}, and efficiently updated (see step 23 in Algorithm 1, and Algorithm 2). The half-space directions are predefined in a constant matrix H ∈ R d×n whose rows are made of d unitary vectors uniformly distributed on an n-dimensional unit hypersphere e.g. as in [START_REF] Marsaglia | Choosing a point from the surface of a sphere[END_REF]. If needed, H can be appended with other specific directions such as those related to the canonical basis, for instance.

Then, Z+ ⊃ Z + is obtained by the iterative update called at step 23 in Algorithm 1. The iteration itself is described in Algorithm 2. For each computed zonotope Z (like Z i in ( 12)), it consists in computing a tight inflation of the bounding vector b so that Z ⊂ Z+ is satisfied. This is basically achieved using the support point and interval hull properties of zonotopes, as expressed in vector form in the body of Algorithm 2.

Algorithm 2 Tight update of half-space bounds of polytope Z+ = {z | Hz ≤ b} to include the zonotope Z = ⟨c, R⟩

1: function b = UpdateBound(H, b, ⟨c, R⟩) 2: b ← max(b, Hc + |HR|1) ▷ 1 is a col. vec. of ones 3: return b
Finally, once the test implemented in Algorithm 1 successfully terminates, the iterative calls of Algorithm 2 result in an intersection of halfspaces, i.e. a polytope, bounding all the trajectories starting from the N -RPI set P. Note that, by definition 2, all these trajectories also get back to P in at most N steps, no matter how the specified attacks and the bounded disturbances change in time.

Robust and Resilient MPC

For the implementation of the robust and resilient control scheme introduced in Section 4, the objective is to find a sufficiently small N -RPI to ensure that the difference between the system dynamic (1) and the nominal model ( 5) meets an acceptable control accuracy (i.e. the required level of precision during normal operation). Moreover, the smaller the N -RPI is, the smaller the trajectory bounding set Z + will be, so leading to less restrictive tightened constraints X = X ⊖ Z + and Ū = U ⊖ KZ + for the MPC control of the nominal model. Also, the larger the N -RPI set is, the larger will be the allowable disturbance bounds and/or the maximum number M of sample attacks as in (3). Thus, for a given scenario, the parameterized (by α) family of sets F P = {αP | α ≥ 1} is considered and the objective is then to find the smallest α such that the set αP is N -RPI for the considered system and scenarios. F P is considered because any linear inflation of an RPI set by a scalar factor greater than one is also an RPI set. As a result, the proposed Resilient and Robust MPC can be set as follows: a) Compute a stable K (e.g. as an LQ solution), b) Compute P as described in [START_REF] Rakovic | Invariant approximations of the minimal robust positively invariant set[END_REF], c) Find α s.t. αP is N-RPI ( §5.1) and compute Z+ ( §5.2), d) Compute X and Ū (polytopic tightened constraints), e) Compute S as a solution of the Lyapunov equation:

(A + BK) ⊤ S(A + BK) -S = -(Q + K ⊤ RK), (13) f) Compute β s.t. the terminal set Xf = {x ∈ R 2 | ||x|| 2
S < β} satisfies Xf ⊂ X and K Xf ⊂ Ū. Note that Algorithm 1 can be called within a dichotomy search to minimize α at step c). Also, e) and f) are focused on the approach concretely used in section 6 to obtain a final weight and terminal set satisfying the requirements of Assumption 3, where the invariance of Xf as in f) is ensured by (13).

NUMERICAL EXAMPLE

Consider a second order system modelling the attack-free dynamic of a CPS as given by (1) and as explained in section 3 with: With those parameters, an ε-approximation of the minimal RPI set P is first obtained using the formula :

A = 1.05 0.5 -0.6 1.1 , B = 1 1 , D = 1 0.5 With the state constraint X = {x ∈ R 2 | ||x|| ∞ ≤ 20}, the control input constraint U = {u ∈ R | ||u|| ∞ ≤ 7},
P(s, γ) = (1 -γ) -1 s-1 i=0 (A + BK) i DW,
as proposed in [START_REF] Rakovic | Invariant approximations of the minimal robust positively invariant set[END_REF], with γ = 0.01 and s = 9. Then, the approach proposed in Section 5 ensures that the set Z = 3P (α = 3) is N -RPI with a maximal number M = 4 of attacks over an N = 10 steps horizon. Algorithm 1 runs in 0.45 seconds on a i5-5200U CPU. To illustrate the use of the obtained N -RPI set, a simulation was done using the attack scenario reported in Fig. 3. To implement the MPC controller, S is computed as a solution of the Lyapunov equation ( 13). This results in S = 2.2567 -0.6841 -0.6841 1.4846 . The terminal set is designed as an ellipsoid X f = {x ∈ R 2 | ||x|| 2 S < β} where β = 0.926 was computed to obtain the biggest set possible while verifying the assumption 3. Set computations were done using the CORA toolbox [START_REF] Althoff | Cora 2021 manual[END_REF]), using zonotope to obtain an RPI and then polytopes for the constraints of the MPC problem. The terminal set was inner-approximated by a polytope for convenience of the MPC implementation.

The online MPC problem ( 8) is formulated as a quadratic programming problem and solved using the function mp-cActiveSetSolver.m provided in the Model Predictive Control Toolbox of Matlab. Starting from an initial state x 0 = [1.4 -3] ⊤ , Fig. 4 shows the tube trajectory of the system along with the initial feasible set. To better observe the resilient set, Fig. 5 shows the trajectory of the difference z k , as well as the sets Z, Z+ and P.

CONCLUSION

A new resilient and robust control scheme is proposed for a class of CPS subject to state and input constraints, unknown but bounded disturbances and possibly repeated time-limited Denial of Service (DoS) attacks cutting an acknowledgement-based communication between the controller and smart 2 actuators. The proposed tube-based MPC scheme was proved to be robust to unknown-butbounded disturbances and jointly resilient to the DoS attacks. To achieve this, a new robust and resilient invariant set, namely µ-RPI set, was introduced in this paper. Finally, an algorithm was proposed to analyze and evaluate the maximal number M of sample DoS attacks over the MPC time horizon N such that the property of µ-step resilience is still robustly ensured. This work paves the way for further investigations. In particular, other common attacks such as false data injection or replay attacks may also be considered. Moreover, the original µ-RPI set introduced in this paper is not specific to an MPC context and may be adapted to other control framework. Special attention may be paid on combining µ-RPI sets with other resilient metrics such as the critical time recently proposed in [START_REF] Perodou | Critical-time analysis of cyber-physical systems subject to actuator attacks and faults[END_REF], all put together to move toward a robust and resilient control of cyber-physical systems.

Fig. 2 .•

 2 Fig. 2. Illustration of Algorithm 1 • Only the scenarios where the first attack happens at the first of N consecutive samples is tested. This is a direct consequence of Assumption 2, where the considered class of attacks satisfies (3) for timehorizon N . • If, for i sample attack instances with the last attack happening at step k a , the system gets back to the RPI set in a total of k b steps, then, to test a scenario with i + 1 attack instances where the first i attacks are the same, testing only the scenarios where the last attack happens between k a + 1 and k b -1 is sufficient.

  and the disturbance set W = {w ∈ R | ||w|| ∞ ≤ 0.1}. For the MPC problem, Q is the identity matrix, R = 0.1, and the prediction horizon is set to N = 10. Then, K is calculated as a solution of the LQ problem with the previous costs.

Fig. 5 .

 5 Fig. 3. DoS attack sequence for 30 times. The blue area denotes the DoS activation times

  b)

	34:	else
	35:	test ← F alse
	36:	end if
	37:

ACKNOWLEDGEMENTS

This study has been carried out with financial support from the French National Research Agency (ANR) in the framework of the Investments for the Future, Programme IdEx Bordeaux-SysNum (ANR-10-IDEX-03-02). The financial support from the "Région Nouvelle-Aquitaine" is also gratefully acknowledged.