C B Faye

T Amodeo

E Fréjafon

N Delepine-Gilon

Christophe Dutouquet
email: christophe.dutouquet@ineris.fr

Sampling considerations when analyzing micrometric-sized particles in a liquid jet using laser induced breakdown spectroscopy

Keywords: suspensions of micrometric-sized particles, LIBS, liquid jet, NOAA (Nano Objects, their aggregates and their agglomerates)

Pollution of water is a matter of concern all over the earth. Particles are known to play an important role in the transportation of pollutants in this medium. In addition, the emergence of new materials such as NOAA (Nano Objects, their Aggregates and their Agglomerates) emphasize the need to develop adapted instruments for their detection. Surveillance of pollutants in particulate form in waste waters in industries involved in nanoparticle manufacturing and processing is a telling example of possible applications of such instrumental development.

The LIBS (Laser-Induced Breakdown Spectroscopy) technique coupled with the liquid jet as sampling mode for suspensions was deemed as a potential candidate for on-line and real time monitoring.

With the final aim in view to obtain the best detection limits, the interaction of nanosecond laser pulses with the liquid jet was examined. The evolution of the volume sampled by laser pulses was estimated as a function of the laser energy applying conditional analysis when analyzing a suspension of micrometric-sized particles of borosilicate glass. An estimation of the sampled depth was made. Along with the estimation of the sampled volume, the evolution of the SNR (Signal to Noise Ratio) as a function of the laser energy was investigated as well. Eventually, the laser energy and the corresponding fluence optimizing both the sampling volume and the SNR were determined. The obtained results highlight intrinsic limitations of the liquid jet sampling mode when using 532 nm nanosecond laser pulses with suspensions.

A C C E P T E D

M A N U S C R I P T

Introduction

Pollution of air and water are major concerns all over the world. At European union level, the importance of the issue of water pollution is stressed by the adoption of a framework directive aiming at achieving a good ecological and chemical status of waters by 2015. Though the water framework directive does not specifically target particles, these are known to play a major role in pollutant transportation. Furthermore, the issue of particle release in liquids could become even more compelling with the advent of new objects such as NOAA (Nano Objects, their Agglomerates and Aggregates). Surveillance of pollutants in particulate form possibly released in waste waters by industries involved in nanoparticle manufacturing and processing is a telling example of applications requiring instrumental developments. The above considerations emphasize the need to develop a technique adapted to particle detection allowing in-situ and real-time elemental identification and mass concentration determination of particles suspended in liquids. The LIBS (Laser-Induced Breakdown Spectroscopy) technique was deemed as a potential candidate to tackle such issue [START_REF] Hahn | Laser-Induced Breakdown Spectroscopy (LIBS), Part I: Review of Basic Diagnostics and Plasma-Particle Interactions: Still-Challenging Issues Within the Analytical Plasma Community[END_REF][START_REF] Hahn | Laser-Induced Breakdown Spectroscopy (LIBS), Part II: Review of Instrumental and Methodological Approaches to Material Analysis and Applications to Different Fields[END_REF].

LIBS consists in focusing a powerful laser pulse on a material (solid, liquid, gas, aerosol, nanoparticle suspensions) whose elemental composition is to be determined. The strong heating of the sample at the focusing spot leads to the ignition of a hot and luminous transient ionized gas called plasma. Plasma light contains the signature of all the chemical elements the interrogated material is made of. This signature is read by sending the emitted light through a spectrometer equipped with a detector. The LIBS signal presents itself as an optical emission spectroscopy spectrum displaying lines corresponding to the detected elements. The features of LIBS (all optical technique, fast, non intrusive, no sample preparation, remote and stand-off analysis possible) do make it a promising analytical chemistry technique intended to be operated at industrial or environmental sites. Current elemental analysis techniques require time consuming procedures involving several steps such as on-site sampling, preparation back to the laboratory prior to the analysis itself.

LIBS has already been applied to particle detection in liquids either for environmental or industrial purposes [START_REF] Ito | Determination of colloïdal iron in water by laser-induced breakdown spectroscopy[END_REF][START_REF] Knopp | Laser-Induced Breakdown Spectroscopy (LIBS) as an analytical tool for the detection of metal ions in aqueous solutions[END_REF][START_REF] Nakamura | Determination of an iron suspension in water by Laser-Induced Breakdown Spectroscopy with two sequential laser pulses[END_REF][START_REF] Haisch | Characterization of colloïdal particles by Laser-Induced Plasma spectroscopy (LIPS)[END_REF][START_REF] Bundschuh | Determination of size concentration and elemental composition of colloids with Laser-Induced Breakdown Detection/Spectroscopy (LIBD/S)[END_REF][START_REF] Yun | Selective determination of Europium (III) oxide and hydroxide colloids in aqueous solution by Laser-Induced Breakdown Spectroscopy[END_REF][START_REF] Hotokezaka | Selective and in-situ determination of carbonate and oxide particles in aqueous solution using Laser-Induced Breakdown Spectroscopy (LIBS) for wearable information equipment[END_REF][START_REF] Yamaguchi | Direct quantitative analysis of particulate aluminum suspended in water using Laser-Induced Breakdown Spectroscopy[END_REF]. Different sampling modes such as cuvette [START_REF] Knopp | Laser-Induced Breakdown Spectroscopy (LIBS) as an analytical tool for the detection of metal ions in aqueous solutions[END_REF][START_REF] Haisch | Characterization of colloïdal particles by Laser-Induced Plasma spectroscopy (LIPS)[END_REF][START_REF] Bundschuh | Determination of size concentration and elemental composition of colloids with Laser-Induced Breakdown Detection/Spectroscopy (LIBD/S)[END_REF][START_REF] Yun | Selective determination of Europium (III) oxide and hydroxide colloids in aqueous solution by Laser-Induced Breakdown Spectroscopy[END_REF][START_REF] Hotokezaka | Selective and in-situ determination of carbonate and oxide particles in aqueous solution using Laser-Induced Breakdown Spectroscopy (LIBS) for wearable information equipment[END_REF][START_REF] Yamaguchi | Direct quantitative analysis of particulate aluminum suspended in water using Laser-Induced Breakdown Spectroscopy[END_REF], liquid jet [START_REF] Ito | Determination of colloïdal iron in water by laser-induced breakdown spectroscopy[END_REF][START_REF] Nakamura | Determination of an iron suspension in water by Laser-Induced Breakdown Spectroscopy with two sequential laser pulses[END_REF][START_REF] Haisch | Characterization of colloïdal particles by Laser-Induced Plasma spectroscopy (LIPS)[END_REF] and deposition on filters (following ultrafiltration) [START_REF] Haisch | Characterization of colloïdal particles by Laser-Induced Plasma spectroscopy (LIPS)[END_REF] were coupled with nanosecond duration laser pulses at 532 nm (for most of them) testing various energies, focusing and collection paths. The analyses were carried out on particles of various chemical natures with sizes ranging from a few hundreds of nanometers to a few tens of micrometers.

Detection limits were found on average of about a few tens to a few hundreds of micrograms per liters.

Continuous monitoring of liquid samples like that needed in the above mentioned possible application is facilitated using the liquid jet configuration. Two recent papers have shed light on peculiarities encountered

A C C E P T E D M A N U S C R I P T

when dealing with particle detection using the liquid jet. Snyder et al. [START_REF] Snyder | The detection of palladium particles in proton exchange membrane fuel-cell water by Laser-Induced Breakdown Spectroscopy[END_REF] have emphasized the need for a careful data processing and rejection of spectra free of signal where no particles are sampled, thereby highlighting the difference between a supposed homogeneous solution and a suspension of particles with their discrete nature.

The upper size limit of complete vaporization was investigated by Diaz Rosado et al. [START_REF] Rosado | Effect of particle size on laser-induced breakdown spectroscopy analysis of alumina suspension in liquids[END_REF] demonstrating incomplete vaporization of micrometric-sized alumina particles above a diameter of around 800 nm.

In this paper, the attention was focused on the volume sampled by the laser pulses and on the SNR (Signal to Noise Ratio) of the LIBS signal. Both sampling volume and SNR evolutions were examined as a function of the laser energy. The obtained sampling volume was compared to the maximum accessible volume calculated considering the liquid jet thickness times the estimated laser beam diameter right where laser pulses strike the liquid jet. A "sampled depth" representing the layer thickness within which the analyte is vaporized thereby contributing to the LIBS signal was estimated. Eventually, a laser energy and an estimation of the corresponding fluence were found that maximize both the sampling volume and the SNR, thereby designating the optimal laser energy / fluence when irradiating a suspension with laser pulses of 5 ns duration at 532 nm. The limitations inherent to the coupling of the liquid jet with the LIBS technique when analyzing micrometric-sized particles are presented and discussed.

Experimental set-up

The liquid sampling mode

The liquid jet sampling mode (designed by IVEA solution French supplier of LIBS systems) has been retained for the experiments on the detection of particles in liquids. The suspension is flowed through a flexible tube connected to a capillary using a peristaltic pump. It is then ejected from the capillary in the form of a stable liquid jet with a 1 mm thickness at a flow rate of 150 mL min -1 . Laser pulses were focused on the jet around 10 mm downstream its output. The IVEA system also includes air sweeping in front of the focusing lens to avoid droplet deposition caused by splashing. Micrometric displacements of the liquid jet was assured using two linear translation stages allowing its precise positioning with respect to the focused laser pulses.

LIBS experimental set-up

The LIBS experimental set-up is presented on figure 1. Laser pulses with energies varying from 3 to 60 mJ and durations of 4 ns fired by a 532 nm Q-switched laser (Quantel Brilliant) operating at 20 Hz were first spatially expanded to three times the initial beam diameter (6 mm). They were then directed through a half wave

A C C E P T E D M A N U S C R I P T
plate and a Glan polarizer for energy attenuation. A dichroïc mirror fully reflective at 532 nm but otherwise transparent directed the beam toward a 100 mm focal length lens allowing focusing laser pulses on the liquid jet.

The mirror was positioned so as to make a 90 degree angle between the optical axis of the beam expander and that of the focusing lens. Such optical path makes light collection along the focusing axis possible, given that the dichroïc mirror is transparent to almost all wavelengths but that of the laser. The LIBS signal was collected along the focusing axis through the 100 mm focal length lens and focused on the entrance of a 1 mm core diameter optical fiber (Ocean Optics) using a 35 mm focal lens. Optical emission spectra were retrieved from a spectrometer (Horiba Jobin Yvon, model iHR320) equipped with a gated intensified CCD camera (Andor iStar model DH734-18F-03). Time resolved measurements are made possible with such detector by triggering signal recording when the laser pulses are fired and recording over the required duration by selecting the appropriate integration time (with accessible time resolution of around 2 ns). Another iCCD camera (same model as above)

equipped with UV lens objective (Nikon model 78 mm F/3.8) was positioned on the axis perpendicular to that of laser focusing to allow monitoring the position of the liquid jet. The focal spot size had to be estimated to evaluate the laser fluence. The estimation was made assuming a perfect Gaussian beam profile and considering the beam section at 1/e 2 of the maximum pulse intensity using marks drilled in a copper plate for different energies as explained in section (4.3). Crater diameters were measured using a Zeiss microscope (Axio Imager M1m).

Foreword

Particle features

Duke scientific (dry glass NIST traceable size standards, 9000 Series Borosilicate Glass) certified particle size standards with lots of particles with nominal diameters of 2, 5, 8 and 10 µm were used in this work.

Four test suspensions containing particles with the aforementioned characteristics and concentrations of around Particle features were eventually found in good agreement with the data sheets provided by the manufacturer even after having been stirred.

Calculation of the sampling volume

As emphasized by Hahn et al. [START_REF] Hahn | Discrete particle detection and metal emissions monitoring using laser-induced breakdown spectroscopy[END_REF] in previous work, the analysis of particles using the LIBS technique requires adapted data processing. Either suspended in the air or in a liquid, particles have a discrete nature. This means that not all the laser shots lead to particle detection. Under these conditions, the usual ensemble-averaging technique aiming at improving the signal-to-noise ratio by systematically adding the signal originating from several laser shots in one final single spectrum is no longer adapted. An alternative method named conditional analysis was proposed by Carranza and Hahn [START_REF] Carranza | Sampling statistics and considerations for single-shot analysis using laserinduced breakdown spectroscopy[END_REF][START_REF] Carranza | Plasma volume considerations for analysis of gaseous and aerosol samples using laser-induced breakdown spectroscopy[END_REF] in order to improve LIBS signal to noise ratio when dealing with particle detection. It consists in recording individual spectra, each of them corresponding to one unique laser shot. The ratio of the number of spectra having sampled one particle (spectra therefore considered as positive) to that of the total amount of laser shots fired [START_REF] Diwakar | Laser-induced breakdown spectroscopy, for analysis of micro-and nanoparticles[END_REF] yields the PSR (or Particle Sampling Rate).

Spectra free of signal are discarded prior to the addition leading to the final spectrum. The latter is obtained by adding the intensities of all the positive spectra and multiplying the resulting sum by the particle sampling rate (PSR). It was shown by Panne and Hahn [START_REF] Panne | Laser-Induced Breakdown Spectroscopy[END_REF] that the sampling volume could be estimated from the particle sampling rate provided that the latter was not higher than a value of around 10 %. With this condition fulfilled, an estimation of the sampling volume can be inferred from equation 1 below: Equation 1with V S the sampling volume, C the number concentration of particles and PSR the particle sampling rate. This approach was applied in order to assess the volume of liquid sampled when analyzing particles using LIBS.

Results and discussion

Line selection and temporal settings

Given the relative abundances of the various elements constituting borosilicate glass particles (see section 3.1) and the efficiency of the detection system as a function of wavelength, the Ca I line at 422.672 nm (A ki = 21.8  10 7 s -1 , E i (eV) = 0.0 and E k (eV) 2.9325118) was selected for these experiments. The evolution of

A C C E P T E D M A N U S C R I P T

the sampling volume and that of the SNR were estimated through the LIBS signal as a function of laser energy.

Accordingly, time delays and gate widths were adjusted after having analyzed LIBS temporal evolution recorded when irradiating a CaSO 4 solution for each laser energy. Such procedure was motivated bearing in mind that the plasma volume tends to enlarge when increasing the energy. Such enlargement results in a much slower cooling speed thereby temporally shifting the time delays at which the recording of the given line of an element is optimized [START_REF] Fisher | Temporal gating for the optimization of laserinduced breakdown spectroscopy detection and analysis of toxic metals[END_REF][START_REF] Amodeo | On-line determination of nanometric and sub-micrometric particle physicochemical characteristics using spectral imaging-aided Laser-Induced Breakdown Spectroscopy coupled with a Scanning Mobility Particle Sizer[END_REF]. The obtained results are presented on figure 4. The optimized time delays were indeed found slightly shifted from one laser energy to another with values of 0.25, 1, 1, 2, 3, 4, 5, 5 µs for laser energies corresponding to 3, 5, 7, 10, 15, 30, 45, 60 mJ. The gate width was fixed to 5 µs.

PSR calculation

The determination of the PSR prior to that of the sampling volume necessitates appropriate data processing. A two step filtering procedure was utilized as suggested by Asgill and Hahn [START_REF] Asgill | Particle size limits for quantitative aerosol analysis using Laser-Induced Breakdown Spectroscopy: Temporal considerations[END_REF]. A first filtering algorithm was designed to select positive spectra. It was operated using a criterion designating a spectrum as positive when the peak of the line of interest minus its background was found higher than a threshold value defined as n times the noise intensity. The latter is calculated on both sides of the line (in its vicinity) of interest using the standard deviation formula. Several values of n were tested when running the algorithm for a given series of spectra. The selection of the threshold level implying the fixing of the n number was set using the following approach. Obviously, low values of n lead to select spectra that do not correspond to real lines. The algorithm is very likely to be misled if the accumulation of noise on pixels where the line of interest is expected to emerge happens to be higher than the defined threshold value. This eventually results in the recording of the so called false positives. Contrariwise, some of the real positives may be accidentally discarded when choosing high values for n. Thus, the n values were chosen so as to minimize the amount of false positives found when visually inspecting the retained spectra.

A second filtering algorithm was implemented to refine the selection. Indeed, poor repeatability of PSR calculation was obtained when only applying the above defined criterion. In addition, lots of the selected spectra (though not clearly identified as false positives) happen to display a very low signal to background ratio. They made the PSR increase without contributing much to the final signal. Thus, the intensity of the Ca line normalized to the maximum recorded intensity was plotted as a function of the percentage of spectra marked as positives (figure 5). Eventually, the PSR repeatability was found to stabilize when suppressing 30% of the

A C C E P T E D M A N U S C R I P T

spectra marked as positives whose intensities were below 10% of the highest recorded intensity. This is tantamount to neglecting minor hits corresponding to partial ablation of the particles. The PSR is therefore slightly underestimated by not taking into account minor hits.

Evolution of the sampling volume as a function of laser energy

A first series of experiments was carried out with the aim in view to investigate the evolution of the sampling volume as a function of laser energy. Three series of spectra of 1000 laser shots each were recorded for several energies ranging from 3 to 60 mJ. Sampling volume variations as a function of laser energy is presented on figure 6. The particle sampling volume V S was inferred from the PSR according to equation 1 as explained above. The obtained curve can actually be divided into two regions. Up to 15 mJ, the sampling volume is seen to grow larger with increasing energy with a steep slope. Contrariwise, only a very slight expansion is to be observed from 15 up to 60 mJ with a much lower slope. Eventually, the particle sampling volume does not prove to linearly expand from 3 to 60 mJ, almost saturating for a laser energy of around 15-20 mJ.

When analyzing silicon targets under similar experimental conditions by LIBS, namely using laser pulses of 5 ns duration at 532 nm, Milan and Laserna [START_REF] Milan | Diagnostics of silicon plasmas produced by visible nanosecond laser ablation[END_REF] observed a curve of similar shape plotting the evolution of the Si line at 288.158 nm as a function of laser fluence. The LIBS signal was found to saturate reaching a certain fluence value. Plasma shielding was concluded to explain such signal saturation. The same is very likely to apply when interrogating the liquid jet. There exists a laser fluence beyond which the sampling volume does not grow anymore owing to plasma shielding, preventing all of the laser pulse to reach the liquid jet.

The above presented results show how the sampling volume evolves with laser energy. Still, both the sampling volume and the SNR are expected to vary not only as a function of the deposited energy, but also through the laser fluence considering the interaction surface between laser pulses and the liquid jet. A first evaluation of this interaction surface can be made estimating the surface of the laser beam at 1/e 2 of the peak intensity, assuming that the laser beam has an almost perfect Gaussian profile.

As no craters or marks may actually be printed on the liquid jet surface, the latter was replaced by a copper plate according to the following procedure. Images of the liquid jet were taken using the iCCD camera equipped with an objective and positioned perpendicularly to the laser beam. Position of the side-viewed surface of the liquid jet facing the laser beam was recorded identifying the corresponding vertical pixel row. The liquid jet was then replaced by a copper plate. Its side-viewed surface exposed to the laser beam was positioned aligning it with the same vertical pixel row as that used to spot the exposed face of the liquid jet. Craters were then drilled in copper for laser energies identical to those used when probing the liquid jet. Crater diameter (considering the inner rim of the crater on the material surface) evolution as a function of laser energy is shown on figure 7. The experimental points were found to fit quite well with a logarithmic equation below writing: Equation 2Assuming that the laser beam used in this study has a perfect round Gaussian shape, the theoretical formula of the surface of the beam section derived from that of the laser fluence writes: Equation 3where S is the crater surface of a particular point in the beam and S 0 is the surface at 1/e 2 . E Laser is the maximum energy at the center of the beam,  0 represents the 1/e 2 radius and  th the melting threshold fluence.

Identifying equation 2 with equation 3 allows determining S 0 whose value was found of about 6500 µm 2 corresponding to a diameter of around 90 µm. The applied laser fluence value was therefore estimated of around 235 J cm -2 with a laser energy of 15 mJ corresponding to saturation onset on figure 6.

Estimation of the sampled depth

The estimated value of the sampling volume was found to tend to saturate for energies of around 15-20 mJ with an estimated laser-jet interaction surface of about 90 µm. Though it may be difficult to determine a precise value of the sampling depth, it may be envisaged to assess its order of magnitude when probing a liquid jet in the conditions described above with nanosecond laser pulses of wavelength at 532 nm. At least two assumptions must be made to make an approximate calculation. First, the geometry of the sampling volume is not known with accuracy and must therefore be approximated by a cylinder. Second, the surface of the base of the cylinder may be approximated by S 0 assuming that the interaction surface between the laser pulses and the liquid jet is close to that value. Eventually, an estimation of the sampled depth h of around 7-10 µm may be drawn from the above calculated values of V s (when saturation occurs, see figure 6) and S 0 . It proves to be much lower than the liquid jet thickness L ( 1 mm) as h / L << 1. This indicates that an in-depth sampling is not achieved under such experimental conditions.

Considering estimations of thermal and optical penetration length compared with the above sampling depth tends to substantiate the hypothesis of plasma shielding. Indeed, depending on the material nature, the area affected by the laser pulses may be dominated either by thermal or optical processes [START_REF] Vladoiu | The dependence of the ablation rate of metals on nanosecond laser fluence and wavelength[END_REF]. The thermal penetration length may be inferred from the following expression: Equation 4where  represents the thermal diffusivity and  the laser pulse duration. Given that the laser pulse duration is of around 4 ns and considering the thermal diffusivity of pure water   1.4  10 -7 m 2 s -1 [START_REF] James | The thermal diffusivity of ice and water between -40 and +60 °C[END_REF], a rough estimation of the thermal penetration length was found of around a few tens of nanometers (L th  25 nm). The optical penetration length is contrariwise much higher than the thermal depth given that water is transparent to the laser wavelength at 532 nm (absorption coefficient   3.5  10 -4 cm -1 [START_REF] Irvine | Infrared optical properties of water and ice spheres[END_REF] and therefore L opt >> 1mm). The liquid jet heating is therefore far more dependent on optical properties than on thermal features. Accordingly, volumic heating prevails over surfacic heating given that L opt >> L th . Still, the sampled penetration length verifies the following inequality when compared with its thermal and optical counterpart: Equation 5This tends to support the hypothesis that the laser pulses can no longer pass through the liquid jet, being blocked by the nascent plasma growing in the liquid.

Comparison of the sampling volume with possible maximum sampling volume

Comparing the measured sampling volume with say, the volume corresponding to the intersection of the laser beam with the liquid jet can help getting a better insight of the efficiency of laser sampling. The latter volume delimitates the largest amount of liquid possibly sampled by the laser (assuming that water sampling is only mostly achieved through laser interaction and not assisted by plasma vaporization). This maximum sampled volume may be approximated multiplying the liquid jet thickness (L  1 mm) by the roughly estimated laser-jet interaction surface S 0 (6500 µm 2 ). The so calculated volume MV S (possible maximum sampling volume) equaling 10 -6 cm 3 proves to be much larger than the sampling volume V S which tends to saturate beyond 15-20 mJ with values of around 6 10 -8 cm 3 (V S / MV S << 1). Thus, the estimations of the sampling volume V S and depth h and their comparison with a possible maximum sampling volume MV S and the liquid jet thickness L respectively make it clear that the area sampled by the laser remains confined near the liquid jet surface. Such

A C C E P T E D M A N U S C R I P T

limitation has important consequences for the analysis of particles using a liquid jet in terms of sampling rate.

Assuming a homogeneous distribution of particles within the liquid jet volume of cylindrical revolution, these results suggest that a large part of the particles will not be analyzed. Indeed, most of those located beyond the estimated sampling depth are not to be taken into account, thereby limiting the sampling rate using such configuration.

LIBS signal evolution as a function of particle sizes

When tackling the issue of the evolution of the LIBS signal-to-noise ratio as a function of various experimental parameters, the dependence of signal response as a function of particle sizes was broached first.

Four suspensions with the same number concentration of 1.47 x 10 6 particles / cm 3 containing particles with monodisperse distributions of 2, 5, 8 and 10 µm were prepared. These were probed by LIBS recording series of individual spectra of 1000 laser shots with energy of 15 mJ. Positive spectra or particle hits were identified following the above described procedure (section 4.2). The same total number of spectra was selected among the positive spectra for each probed suspensions, and therefore each particle size. Four final spectra were obtained by ensemble averaging an equal number of particle hits. The dependence of peak to base ratio as a function of particle volume is presented on figure 8. The obtained graph shows that the LIBS signal does not vary as a function of particle size. As the LIBS signal is known to be proportional to the ablated mass, a linear dependence would be expected when carrying out these experiments. Two possible explanations may be suggested. Either such effect is the consequence of self-absorption or of non-complete vaporization of the particles, like that observed by Carranza and Hahn [START_REF] Carranza | Assessment of the upper particle size limit for quantitative analysis of aerosols using Laser-Induced Breakdown Spectroscopy[END_REF] when analyzing particles in air using LIBS. A similar result has recently been found by Diaz Rosado et al. [START_REF] Rosado | Effect of particle size on laser-induced breakdown spectroscopy analysis of alumina suspension in liquids[END_REF] when analyzing micrometric-sized alumina particles in liquids using LIBS. A careful examination of the recorded data showed that non complete vaporization was very likely to be the cause of signal saturation. The graph of figure 8 shows likewise that borosilicate particles were not entirely vaporized either for sizes ranging from 2 to 10 µm.

SNR evolution as a function of laser energy

As explained above, the estimation of the sampling volume is of primary importance to pinpoint the intrinsic limitations of the coupling of LIBS with the liquid jet as a sampling technique for particle analysis and therefore the detection limits inherent to such approach. Investigations on the SNR must be carried out along with those on the sampling volume in order to assess the LIBS response. The evolution of the SNR as a function

A C C E P T E D M A N U S C R I P T

of laser energy (see figure 9) was extracted from the very spectra recorded when analyzing the sampling volume (section 4.2). As already stated when examining the graph of figure 6 the curve depicting the evolution of the SNR as a function of laser energy may be divided into two regions. The SNR linearly rises from 3 mJ up to 15-20 mJ. Such rise is likely to be accounted for by the expansion of the sampling volume as depicted on figure 6.

All like the sampling volume, the SNR happens to saturate for energies beyond 15-20 mJ for the same reason than that mentioned when scrutinizing figure 6, namely plasma shielding.

Optimizing the sampling volume and the SNR through proper laser fluence settings

The optimization of the LIBS signal (under the experimental conditions described in this paper) implies that of the sampling volume along with that of the SNR. Both the SNR and the sampling volume were found to saturate beyond a laser energy of around 15-20 mJ corresponding to an estimated fluence of around 235 J/cm 2 .

This result shows that systematically applying high laser fluences is no use, as neither the sampling volume nor the SNR increase with increasing energy. Increasing energy will only enhance mechanical ablation (caused by the action -reaction principle, the expanding plasma pushing back the liquid jet) within the liquid jet thereby making the inevitable problem of splashing gaining in importance.

Conclusion & perspectives

LIBS is deemed as a possible tool for in-situ analysis of particles for applications such as wastewater monitoring. With the aim in view to evaluate its potentialities in this context, suspensions of monodisperse borosilicate glass particles with sizes of 2, 5, 8 and 10 µm were analyzed coupling a liquid jet and LIBS.

Nanosecond laser pulses at 532 nm with energies ranging from 3 to 60 mJ were utilized. The Ca I line at 422.672 nm was selected for spectrum recordings, the latter element being one of the major components of borosilicate particles. Thenceforth, the evolution of the sampling volume and of the SNR as a function of the laser energy was investigated. The dependence of the LIBS signal with different particle diameters was examined as well.

The sampling volume (along with the SNR) has been assessed using a conditional analysis procedure close to that implemented by Asgill and Hahn [START_REF] Asgill | Particle size limits for quantitative aerosol analysis using Laser-Induced Breakdown Spectroscopy: Temporal considerations[END_REF] when dealing with the analysis of micrometric-sized particles in air. It was found to linearly rise with laser energy up to 15 mJ whereas it exhibited saturation beyond this value. Such behavior was attributed to plasma shielding. Assuming an interaction surface equaling that at 1/e 2 of the laser peak intensity experimentally estimated drilling marks in a copper plate, an estimation of the sampled depth was made. Along with these, a maximum accessible volume for sampling delimitated by the interaction of

A C C E P T E D M A N U S C R I P T

the laser beam with the liquid jet was calculated. The comparison of the sampled volume and depth with the maximum accessible volume and the liquid jet thickness respectively clearly showed the surfacic nature of the interaction of the laser pulses with the liquid jet. Not all of the liquid jet thickness happens to be probed by the laser beam, thereby suggesting that particles positioned beyond the sampling depth are not taken into account for the analysis. The SNR was demonstrated to display the same tendency as that depicted by the sampling volume with signal saturation beyond 15 mJ. Last but not least, the LIBS signal (peak to base ratio) was found constant as a function of particle diameter, suggesting possible non-complete vaporization of particles. The latter explanation was confirmed by Rozado et al. [START_REF] Rosado | Effect of particle size on laser-induced breakdown spectroscopy analysis of alumina suspension in liquids[END_REF] in an article dedicated to size effects occurring when analyzing particles by LIBS with a liquid jet.

Eventually, all these results point out the precautions of use and limitations of the liquid jet applied to particle analysis with LIBS using nanosecond pulse duration lasers at 532 nm. Indeed, laser fluence is to be set with care given that both the sampling volume and the SNR were seen to saturate. Increasing the laser fluence would only lead to increase splashing without improving the LIBS signal, only favoring mechanical ablation over liquid vaporization. Moreover, quantitative analysis of micrometric-sized particles may be seriously limited by non-complete vaporization of particles. The conclusions drawn from these experiments suggest the use of a liquid sampling technique more adapted to the issue of micrometric-sized particle detection than the liquid jet.

A C C E P T E D M A N U S C R I P T

A C C E P T E D M A N U S C R I P T

10 - 4 g

 4 / cm 3 were prepared using ultrapure water (Millipore water purification system Milli-Q Advantage A10) without surfactant. Particle morphology, elemental composition and sizes were controlled using SEM (Scanning Electron Microscopy, FEI Company, model Quanta 400) analysis. Size distributions of the suspensions were obtained using LDS (Laser Diffraction Spectrometry, Malvern Instruments model Mastersizer 2000). It should be noted that these measurements were performed while continuously agitating the suspensions with a magnetic stirrer, the latter being utilized in the course of the LIBS experiments. SEM images of 5µm sized particles are presented on figure 2. Their nominal diameter was found of around 5  0.5 µm with a spherical shape. Their elemental chemical composition was analyzed through EDX (Energy Dispersive X-ray) showing that Si, Ca and Al are elements found with the largest abundance. Size distributions of the four lots are displayed on figure3.

Figure 1 .

 1 Figure 1. LIBS experimental setup.

Figure 2 .

 2 Figure 2. SEM image of 5 µm-sized borosilicate glass particles.

Figure 3 .

 3 Figure 3. Size distributions of borosilicate glass particles with sizes of 2, 5, 8 and 10 µm obtained using Laser Diffraction Spectrometry (LDS).

Figure 4 .

 4 Figure 4. Temporal evolution of the signal to noise ratio (SNR) corresponding to the Ca I line at 422.672 nm obtained when analyzing a CaSO 4 solution for several energies ranging from 3 to 60 mJ.

Figure 5 .

 5 Figure 5. Distribution of the percentage of spectra as a function of signal intensity normalized to the highest peak obtained when analyzing a suspension of 5 µm-sized particles using an energy of 15 mJ and time delay and gate width of 3 and 5 µs respectively. The vertical line represents the cut-off value under which peaks are no longer taken into account in the PSR calculation.

Figure 6 .

 6 Figure 6. Evolution of the sampling volume as a function of the laser energy recorded for three series of spectra of 1000 laser shots each.

Figure 7 .

 7 Figure 7. Evolution of the diameters of craters drilled on copper as a function of laser energy. The curve has been fitted by a logarithmic expression to deduce the laser beam section at 1/e 2 of the peak intensity.

Figure 8 .

 8 Figure 8. Evolution of the peak to base (P/B) ratio of the Ca line at 422.672 nm as a function of the particle volume for particulates with sizes of 2, 5, 8 and 10 µm recorded obtained when using a laser energy of 15 mJ.

Figure 9 .

 9 Figure 9. Evolution of the signal to noise ratio (SNR) as a function of laser energy obtained using the same set of data as that utilized for figure 6.

Figure 2 AFigure 3 Figure 4 A 20 Figure 5 Figure 6 A 21 Figure 7 Figure 8 A

 23420562178 Figure 2

Figure 9 A 23 HIGHLIGHTS

 923 Figure 9

ACCEPTED MANUSCRIPT

Acknowledgements

This research was supported by the French ministry of ecology, sustainable development and energy. We would like to thank M. Patrice Delalain for his assistance when using SEM and LDS instruments.

A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T


Ana Paula Estrada-Vargas

Ernesto López-Mellado

Jean-Jacques Lesage
email: jean-jacques.lesage@lurpa.ens-cachan.fr

Identification of Partially Observable Discrete Event Manufacturing Systems

Keywords: Black-box identification, Automated manufacturing systems, Interpreted Petri nets, Internal behaviour model. I

This paper deals with black-box identification of discrete event manufacturing systems that are automated using a programmable logic controller (PLC). The behavior of the system is observed during its operation and is represented by a single long sequence w of observed input/output (I/O) signals vectors. The identification method, conceived for addressing large and complex industrial systems, consists of two complementary stages; the first one obtains, from w, the observable part of an interpreted Petri net (PN) model composed of observable places and transitions describing the reactive behavior of the system. Afterwards, w is transformed into a sequence S of transition firings from which a PN model that reproduces S is inferred. This paper focuses on the second stage of the method in which a PN is built by adding non-labeled places and arcs that represent the non-observed behavior of the whole system by assuring the reproduction of w; this technique is based on discovering the causal and concurrent relationships between transitions in S.

INTRODUCTION

Identification methods of discrete event systems (DES) allow building systematically a mathematical model (Petri nets, automata) that describes the behavior of an unknown or ill-known system based on the observation of its evolution. Observations consist of data revealing the system activity: sequences of operations, events, messages, signals etc., and the models allow reproducing the observed behavior.

DES identification has been first addressed as a problem of grammatical inference [START_REF] Gold | Language Identification in the Limit[END_REF] [START_REF] Angluin | Queries and Concept Learning[END_REF] for obtaining finite automata (FA) that represent a given language. Afterwards, Petri net (PN) models have been proposed for coping with more complex systems exhibiting concurrent behavior [START_REF] Hiraishi | Construction of Safe Petri Nets by Presenting Firing Sequences[END_REF].

Three main approaches more specifically conceived for identifying discrete event manufacturing systems have been proposed in recent literature.

The incremental synthesis approach, proposed in [START_REF] Meda-Campaña | A passive method for on-line identification of discrete event systems[END_REF] [5], deals with unknown partially measurable DES exhibiting cyclic behavior. Several algorithms for building interpreted PN (IPN) have been proposed allowing the on-line identification of concurrent DES from output sequences. Although the techniques are efficient, the obtained models may represent more behaviors than those observed.

Other recent method [START_REF] Klein | Fault detection of Discrete Event Systems using an identification approach[END_REF] allows building efficiently a non-deterministic FA (NFA) from a set of input/output sequences, measured from the DES to be identified. The obtained NFA generates exactly the same input/output (I/O) sequences of given length than the observed ones. The method was conceived for fault detection in a model-based approach [START_REF] Roth | An FDI Method for Manufacturing Systems Based on an Identified Model[END_REF] and extended for obtaining an optimal partitioning of concurrent subsystems for distributed fault detection [START_REF] Roth | Black-box identification of discrete event systems with optimal partitioning of concurrent subsystems[END_REF]. The off-line approach based on integer linear programming (ILP) yield free-labeled PN models representing exactly the observed behavior expressed as sequences of events [START_REF] Giua | Identification of free-labeled Petri nets via integer programming[END_REF]. The method is able to handle few short sequences, due to the inherent limitations of ILP regarding its computational complexity. This approach is being explored for other PN classes [START_REF] Cabasino | Identification of Petri Nets from Knowledge of Their Language[END_REF] [START_REF] Dotoli | Real time identification of discrete event systems using Petri nets[END_REF]. A recent stochastic approach allows obtaining timed PN models [START_REF] El Mehdi | Design and identification of Stochastic and Deterministic-Stochastic Petri Nets[END_REF]. Other related works can be found in surveys on identification methods in [START_REF] Estrada-Vargas | A Comparative Analysis of Recent Identification Approaches for Discrete-Event Systems[END_REF] and [START_REF] Cabasino | Model identification and synthesis of discrete-event systems[END_REF]. Furthermore, recent publications on process mining techniques, more suitable for event driven organizational systems than for industrial ones, can be found in [START_REF] Van Der Aalst | Process Mining: Discovery, Conformance and Enhancement of Business Processes[END_REF].

In our approach the problem of identifying partially observable discrete manufacturing systems composed by a controller (a Programmable Logic Controller: PLC) and a plant operating in closed loop is addressed. Both controller's inputs and outputs are sampled for building a single sequence of I/O vectors, which is processed yielding an IPN model.

The aim is to discover, from this observation, how operations of the system are interrelated and construct a concise model which can explicitly show the discovered behavior, in particular, concurrency, synchronization, resource sharing, etc.

Identification of systems in operation involves two important aspects to consider: the system operation and the observation process. Technological issues of both aspects must be considered in the proposed algorithms to construct suitable abstractions. In this paper these issues are addressed by analyzing the observed sequence in order to establish a clear relationship between inputs and outputs of the controller. The proposed method allows building a compact and expressive IPN that is ordinary and safe. It consists of two complementary stages; the first one obtains, from the I/O sequence w, the observable subnet composed by places and transitions labeled with output and input functions respectively [START_REF] Estrada-Vargas | Identification of Industrial Automation Systems: Building Compact and Expressive Petri Net Models from Observable Behavior[END_REF]; during the construction of the model a transition sequence S, which reproduces w, is built.

The paper focuses on the second stage, which allows building efficiently from S, the non-observable part of the model including places (and arcs) ensuring the reproduction of w. The remainder of the paper includes a brief recall of the first stage of the method, and develops the proposed approach for building the non-observable part of the identified PN.

II. BACKGROUND

This section presents the basic concepts and notation of PN and IPN used in this paper. The incidence matrix of

G is C = C + -C -, where C -= [c ij - ]; c ij -= I(p i , t j ); and C + = [c ij + ]; c ij + = O(p i , t j
) are the pre-incidence and post-incidence matrices respectively.

A marking function M : P→ Z + represents the number of tokens residing inside each place; it is usually expressed as an |P|-entry vector. Z + is the set of nonnegative integers.

Definition 2: A Petri Net system or Petri Net (PN) is the pair N = (G,M 0 ), where G is a PN structure and M 0 is an initial marking.

In a PN system, a transition t j is enabled at marking M k if ∀p i ∈ P, M k (p i ) ≥ I(p i , t j ); an enabled transition t j can be fired reaching a new marking M k+1 . This behavior is represented as M k ⎯→ ⎯ j t M k+1 . The new marking can be computed as M k+1 = M k + Cu k , where u k (i) = 0, i≠j, u k (j) = 1; this equation is called the PN state equation. The reachability set, denoted by R(G,M 0 ), of a PN contains all possible reachable markings from M 0 firing only enabled transitions. A PN is said to be 1bounded or safe when

∀M k ∈R(G,M 0 ), ∀p i ∈P, M k (p i )≤1.
Now it is defined IPN [START_REF] David | Petri Nets for Modeling of Dynamic Systems-A Survey[END_REF], an extension to PN that allows associating input and output signals to PN models.

Definition 3 : An IPN (Q, M 0 ) is a net structure Q = (G, V, Σ, Φ, λ, ϕ) with an initial marking M 0 where: G is a PN structure, V = {v 1 , v 2 , ..., v r } is the set of variables, Σ = {α 1 , α 2 ,
..., α s } is the set of events, and Φ = {φ 1 , φ 2 ,..., φ q } is the output alphabet. λ : T→ C × E is a labeling function of transitions, where C={C 1 , C 2 ,…} is the set of variable conditions and E={E 1 , E 2 ,…} is the set of events. In an IPN, a transition t j can be fired if t j is enabled, and if condition C(T j ) is true, when the event in E(T j ) occurs.

ϕ : R(Q,M 0 )→(Z + ) q is an output function, that associates to each marking in R(Q,M 0 ) a q-entry output vector; q=|Φ| is the number of outputs. ϕ is represented by a q×|P| matrix, such that if the output symbol φ i is present (turned on) every time that M(p j ) ≥ 1, thenϕ (i, j) = 1, otherwise ϕ(i, j) = 0.

The state equation is completed with the marking projection Y k = ϕM k , where Y k ∈ (Z + ) q is the k-th output vector of the IPN. Definition 4: A place p i ∈P is said to be observable if the ith column vector of ϕ is not null, i.e. ϕ(•,i) ≠ 0. Otherwise it is non-observable. P = P o ∪ P u where P o is the set of observable places and P u is the set of non-observable places.

III. INPUT-OUTPUT IDENTIFICATION

The problem statement and the main features of the identification method we propose are briefly described. A more detailed presentation can be found in [START_REF] Estrada-Vargas | Identification of Industrial Automation Systems: Building Compact and Expressive Petri Net Models from Observable Behavior[END_REF].

A. Identification of automated DES

In this work we consider systems composed by a Controller (a PLC) and a Plant, denoted as {PLC + Plant}, working in a closed loop. The input signals of the PLC (outputs of the Plant) are generated by the sensors of the Plant. The output signals of the PLC (inputs of the Plant) control the actuators of the Plant.

The identification is made from the point of view of the PLC (Fig. 1). A PLC cyclically performs three main steps: a) Input reading, where signals are read from the sensors; b) Program execution, to determine the new outputs values for the actuators; and c) Output writing, where the control signals to the actuators are set. At each end of the Program execution step, the current value of all Inputs and Outputs (called I/O vector) is captured and, if it is different from the previous one, recorded in a data base. Besides the number of inputs and outputs, the only available data for the identification procedure is a single I/O vector sequence, in which two consecutive vectors are different, whose length depends on the observation duration:

,... ) 3 ( ) 3 ( , ) 2 ( ) 2 ( , ) 1 ( ) 1 ( ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = O I O I O I w
where I(j) and O(j) are respectively the values of the r inputs and q outputs at the j-th PLC cycle.

In order to analyze signals evolution, we compute event vectors, i.e., the difference between two consecutive I/O vectors: E(k) = w(k + 1)w(k) ≠ 0. Each event vector can be decomposed into input and output event vectors:

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = ) ( ) ( ) ( k OE k IE k E in

B. Identification of Partially observable DES

All of the above situations must be taken into account to represent the system dynamics in the model. Besides representing the system's behavior from the I/O vector sequence as an IPN, our purpose in this research is also to provide an identified model as expressive and compact as possible, allowing to represent causal relationships and concurrency of the involved operations. For this purpose we have conceived a two-step method based on the analysis of the relationships between inputs and outputs along the observed behavior represented by the I/O-sequence w. The two steps of the method are the following:

Step1. Discovering the input-output reactive behavior. The observable part of the IPN is built, consisting of sub-graphs named "IPN fragments", composed by observable places and transitions labeled with algebraic expressions onto the set of input variables. Figure 2 shows an IPN model including four fragments.

Step2. Inferring the non-observable part of the IPN. All the observable behaviors being captured in IPN fragments previously built, the sequence w is transformed into a firing sequence S of observable transitions. Non-observable places and arcs, representing an abstraction of the internal and nonobservable behavior of the identified system, are added in such a manner that S (thus w) can be reproduced in the final IPN. In Fig. 2 the non observable part is depicted with dashed lines.

Fig. 2 IPN model including four fragments and non observable nodes

Notice that IPN fragments can be connected, through nonobservable places, only sequentially or in parallel. That's why in section III we will present a technique for discovering causal and concurrency relationships between transitions in S.

The outcome of the first stage of the method is now briefly recalled with the help of a small example in next section.

C. Example 1

The purpose of this example is to show the outputs of the step 1 of the method when a sequence w is processed; it is fully developed in [START_REF] Estrada-Vargas | Identification of Industrial Automation Systems: Building Compact and Expressive Petri Net Models from Observable Behavior[END_REF]. Consider the manufacturing system shown in figure 3 whose function is to sort parcels according to their size. The PLC has 9 inputs, that are signals generated by the sensors of the plant for detecting positions of cylinders (a 0 , a 1 , a 2 , b 0 , b 1 , c 0 , c 1 ) and presence of parcels (k 1 , k 2 ), and 4 outputs, that are signals controlling the actuators of the plant (cylinders A+, A-, B, C). The vector entries correspond to the following distribution: Conveyor 2 (small parcels)

[A+ A-B C k 1 k 2 a 0 a 1 a 2 b 0 b 1 c 0 c 1 ] T .
0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 1 0 1 0 0 1 0 1 0 1 0 0 0 1 0 1 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 1 0 1 0 0 1 1 0 0 0 0 0 0 1 0 1 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 1 0 1 0 0 1 0 1 0 1 0 0 0 1 0 1 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 1 1 0 0 0 0 1 0 1 0 1 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 …
Conveyor3 (large parcels) A+ A- a 0 a 1 a 2 k 2 k 1 B C c 0 b 0 b 1 c 1 1 t 5 t 2 t 6 t - A 4 t B C 3 t + A 7 t ) 1 _ 1 ( ) 1 ( ) ( ) 1 _ 2 ( ) 1 ( ) ( ) ( ) 0 0 0 2 ( ) ( ) 1 _ 0 ( ) 1 ( ) ( ) 1 _ 1 ( ) 1 ( ) ( ) 1 _ 1 ( ) 1 ( ) ( ) ( ) 0 0 0 1 ( ) ( 7 6 5 4 3 2 1 c t a t c b a k t a t b t a t c b a k t • = = • = = • ∧ ∧ ∧ = • = = • = = • = = • ∧ ∧ ∧ = λ λ ε λ λ λ λ ε λ S =
t1 t2 t3 t4 t1 t2 t3 t4 t1 t2 t3 t4 t5 t6 t7 t4 t1 t2 t3 t4 t5 t6 t7 t4 t1 t2 t3 t4 t1 t2 t3 t4 t5 t6 t7 t4 t5 t6 t7 t4 t1 t2 t3 t4 t5 t6 t7 t4 t5 t6 t7 t4 t5 t6 t7 t4 t5 t6 t7 t4 t5 t6 t7 t4 t1 t2 t4 t3 t1 t2 t3 t4 t1 t2 t3 t4 t1 t2 t4 t3 t1 t2 t3 t4 t1

IV. IDENTIFICATION OF THE NON-OBSERVABLE BEHAVIOR

A. Problem (re)statement

The second stage of the method determines the non-observable part of the model consisting of pertinent unlabeled places (and arcs) that rely the fragments and assures that the sequence S (consequently w, which is the actual observed behavior) can be reproduced.

This problem can be stated as follows: given an observable IPN model whose structure is (P obs , T, Pre obs , Post obs ) and a transitions sequence S = t 1 t 2 … t j … ∈ T* reproducing the I/O sequence w, a PN structure (P nobs , T, Pre nobs , Post nobs ) that reproduces S and an initial marking M 0 enabling S must be found. The new PN structure is G=(P, T, I, O) with P= P obs ∪ P nobs , Pre= Pre obs ∪ Pre nobs , Post= Post obs ∪ Post nobs . The PN must be ordinary and safe.

The method proposed herein extracts, from S, precedence and concurrency relations among transitions, which will determine univocally the final structure of the identified model. First, some properties derived from the sequence S are introduced. Afterwards, based on such properties, a technique allowing determining sequential and concurrency relationships among the transitions in S is proposed. Then, the rules for building a net structure rendering the sequential and concurrency relationships are presented.

B. Behavioural properties

Since our construction method is based on the discovering of sequential and concurrency relationships into S, some notions must be defined before introducing the construction procedure of the non-observable behavior.

Definition 5. The relationship between transitions in S that are observed consecutively is expressed in a relation Seq ⊆ T × T which is defined as Seq ={(t j , t j+1 )| 1 ≤ j < |S| }. If (t a , t b ) ∈ Seq, this is denoted by t a <t b .

In a PN model every pair in Seq may in fact be represented differently. If t a <t b , this behavior could be issued from one of two situations in N described in the following definition. Definition 6. Every couple of consecutive transitions t a , t b in Seq can be classified in one of the following situations: -Causal relationship. If the occurrence of t a enables t b . In a PN structure, this implies that there must be at least one place from t a to t b . -Concurrent relationship. If both t a and t b are simultaneously enabled, but t a occurs first and its firing does not disable t b . In a PN structure, this implies that it is impossible the existence of a place from t a to t b . In this case, t a and t b are said to be concurrent, denoted as t a ||t b .

The following notion is the systematic precedence of a transition t j with respect to another transition t k ; it establishes a necessary condition for t j to occur repeatedly. Definition 7. A transition t j is preceded systematically by t k , denoted as t k ∠t j iff t k is always observed between two apparitions of t j in S. By convention, we say that t j ∠t j if t j was observed at least twice in S. The Systematic Precedence Set of t j is a the function SP: T→2 T , that indicates which transitions must be fired to re-enable the firing of t j , i.e. SP(t j )={t k |t k ∠t j }. If t j was observed only once in S, then SP(t j ) = ∅.

The following notion determines straightforward a particular structure from S.

C. Causal and concurrency relationships

Based on the previous definitions we can now determine, for every pair of consecutively observed transitions in S, if the causality or concurrence property is verified. We will therefore deduce some structural properties regarding N. Definition 9. A PN circuit is a path starting and ending in the same node. A circuit is said to be simple if it does not use the same transition more than once, and elementary if it does not use the same place more than once.

Relations between transitions in S can be determined. Proofs are omitted for the sake of brevity; they can be found in [START_REF] Estrada-Vargas | Black-box Identification of Automated Discrete Event Systems[END_REF]. If a couple of transitions (t a ,t b ) in Seq belongs also to CausalR, there must be a place from t a to t b in order to preserve the observed firing order. For the rest of the transition pairs in Seq, we must decide if a place should exist to relate them. Below, we will discuss some cases where the existence of a place can be discarded.

C.1Causal relationship

C.2 Concurrency relationship

If two transitions t a and t b are concurrent, there must not exist a place neither from t a to t b nor from t b to t a ; otherwise, the firing of one would constrain the firing of the other one.

Definition 11. The set of all pairs of concurrent transitions is called

ConcR={(t a ,t b )| t a ||t b }.
If the sequence w is complete, (consequently, S) i.e., if it shows all of the possible behavior of the observed system, we can find concurrency between transitions that are not in a causal relation, as showed in the next proposition. In the sequence of Example 2, t 3 and t 4 fulfill the conditions to be concurrent. They are the only concurrent transitions ConcR={(t 3 ,t 4 )}.

It is well known that, in practice, the sequence w is not complete since it is not possible to assure that the whole behavior of a system has been observed in a finite time. The condition of Proposition 5 is therefore very restrictive, since it requires the observation of all possible behaviors; it could lead to the construction of incorrect models in case of incomplete sequences. Then, less constraining rules to find concurrence must be considered. Next, we present several properties which allow us to identify pairs of transitions which must be concurrent in the identified net N.

First, the notion of Sequential Independence, which is a characteristic of concurrent transitions, is introduced. Later, the propositions to find concurrency will be presented. In the example, we know that t 3 ||t 4 , however, condition b) of Proposition 5 is not fulfilled because SP(t 4 ) = {t 4 }. The SP(t j ) is a singleton, meaning that t j belongs to more than one elementary circuits and then Proposition 5 does not allow to find concurrent transitions to t j . But if t j is included in the SP of other transitions, we can find some concurrency relations, as shown in the next proposition. Proposition 6. Let t a and t b be two transitions in S that have been observed consecutively in both orders (t a < t b and t b < t a ). t a ||t b if t a and t b a) are Sequentially Independent, and b) there exists a transition t k such that t a ∠t k (t a ∈SP(t k ))

and t b ∠t k (t b ∈SP(t k )) In the current example t3 and t4 fulfill the conditions of this proposition, i.e. t 3 <t 4, t 4 <t 3, both transitions are sequentially independent and belong to SP(t 2 ).

If concurrent transitions do not belong to synchronized threads, conditions of the next propositions help us to find a subset of concurrent transitions in which their firings do not depend from the occurrence of another transition t k .

Proposition

) SP -1 (t a ) ≠ ∅, ∀t j ∈SP -1 (t a ), t j || t b .

D. Building the non-observable PN

The computed causal and concurrency relations from sequence S are used to infer internal evolutions of the system by computing non-observable places of the net. Definition 13. The set Seq'= (Seq \ CausalR) \ ConcR contains the set of transition pairs (t a ,t b ) which have been observed consecutively, but are not in a causal relationship or in a concurrency relationship. If Seq' ≠ ∅, there are two possibilities for the remaining transition pairs (t a ,t b ) in Seq': a) They are input and output transitions respectively of a place with several input and output transitions. b) They are concurrent, but w (thus, S) is not complete enough to find such a relationship. Since our goal is to approximate as much as possible the language generated by the net N to the observed sequence S, we assume that if we have observed two transitions consecutively (t a <t b ) but none of the previous propositions have determined that they are concurrent, thus the firing of t a has enabled t b . This is made in order to preserve in N the firing order observed in S. Then, a place will be added from t a to t b ; this is denoted by [ For the parcels sorting system of Fig. 1, the PN structure as well as the computed initial marking (which is computed by allocating tokens enabling the sequence S) is shown in Fig. 5.a. Observe that in Fig. 4 observable places [t 2 , t 3 ] and [t 6 , t 7 ] already exist. By adding non-observable places to such a model and deleting implicit places [t 2 , t 3 ] and [t 6 , t 7 ], we obtain the IPN model shown in Fig. 5.b, which reproduces w.

E. Token flow verification

The sequence w may not have shown enough combinations allowing determining concurrency. If w were complete, all the concurrent and sequential behavior could be found and represented, according to Proposition 5. However, since we know that w could not be complete, in order to approximate the language of N to S as much as we can, we have considered that if two transitions have not been declared as concurrent, they must be in a sequential relationship. However, if the transitions are actually concurrent, the sequential consideration could lead us to create arcs or places in the model constraining too much the behavior of the system and do not allow the firing of S. Now, we present some notions to verify if added places do not interfere in the correct reproduction of S.

F. Summary of the method

All the algorithms described in this section to construct the non-observable part of the IPN are summarized in the following polynomial-time procedure.

Algorithm 1. Non-observable behavior construction

0 0 0 1 : 1 c b a k t ∧ ∧ ∧ 0 0 0 2 : 5 c b a k t ∧ ∧ ∧ 1 : 2 a t ↑ 2 : 6 a t ↑ - A 0 : 4 a t ↑ B C 1 : 3 b t ↑ + A 0 : 7 c t ↑ b)

G. Implementation issues

The algorithm has been implemented and tested on many examples exhibiting diverse situations using the following scheme: a PN is built with the help of a PN editor (PIPE), and then S is created by firing (enabled) transitions randomly. After the processing of S a PN model is identified, and then coded in XML, to be displayed again with PIPE. As an example, an interesting feature of the method is showed. The identification method applied to the following sequence: The complete method including both stages has been implemented as a software tool that process a sequence w and yields an IPN, reproducing exactly w. It has been tested on sequences obtained from real manufacturing systems. For a detailed description please consult chapter 5 in [START_REF] Estrada-Vargas | Black-box Identification of Automated Discrete Event Systems[END_REF].

S = t

ACKNOWLEDGEMENT

The first author has been supported by CONACYT, Mexico, Grant No. 50312 and Région Ile de France.

CONCLUSIONS

A novel approach for identifying partially observable DES systems has been proposed. This method is composed of two steps: the first one is devoted to build the observable reactive behavior [START_REF] Estrada-Vargas | Identification of Industrial Automation Systems: Building Compact and Expressive Petri Net Models from Observable Behavior[END_REF]; the second step, devoted to infer the nonobservable behavior, has been described in this paper.

Although in literature there are techniques addressing the reformulated problem dealing with sequences of transitions, in this paper neither the cycles nor the language generated by the system are known a priori, and this method can handle very long sequences efficiently.

The complete method copes with complex automated DES because it takes into account technological characteristics of both actual controlled systems and the observation process. Implemented as a software tool, the method has been tested on experimental case studies which are very close to actual industrial discrete event processes. The performed tests reveal the efficiency of the methods when sequences formed by thousands of input-output vectors are processed in few seconds [START_REF] Estrada-Vargas | Black-box Identification of Automated Discrete Event Systems[END_REF].

Since the approach is black-box, the obtained models represent the observed behavior; consequently, when the observation is made for a long time, the IPN approximates closely the actual behavior; afterwards the model can be completed using available knowledge on the process.

Current research addresses the problem of discovering indirect causal relationships of transitions that have not been observed consecutively by determining the t-invariants from S.

Definition 1 :

 1 An ordinary Petri Net structure G is a bipartite digraph represented by the 4-tuple G = (P, T, I, O) where: P = {p 1 , p 2 , ..., p |P| } and T = {t 1 , t 2 , ..., t |T| } are finite sets of vertices named places and transitions respectively; I(O) : P × T → {0,1} is a function representing the arcs going from places to transitions (from transitions to places).

Fig. 1 .

 1 Fig. 1. {PLC + Plant} compound and identification procedure.

  which there are three possible observable situations: (j) ≠ 0 and OE(j) ≠ 0: An input change has provoked directly an output change; consequently, a state evolution. The I/O causality is observed at the same PLC cycle. b) IE(j) ≠ 0 and OE(j) = 0: Only the input change is observed in the PLC cycle. It could mean that either i) The output change will appear in the next cycle, or ii) The controller is not sensitive to the input change, or iii) It provoked a change of state that is not observable through a change of output. c) IE(j) = 0 and OE(j) ≠ 0: Only the input change is observed in the PLC cycle. It could be a consequence of the previous situation (case i) or a programmed behavior considering only input levels or timers from the state reached at the event j-1.

Figure 3 .

 3 Figure 3. System's layout A sequence of 216 observed I/O vectors has been processed (only the twenty first vectors are shown below). The IPN fragment, the transition labeling functions, and the transition sequence S shown in figure 4 are obtained. Notice that in this example only one fragment is found.

Figure 4 .

 4 Figure 4. IPN fragment, transition labeling, and sequence of transitions resulting from the first stage of the identification method.

Fig. 5 .

 5 Fig. 5. a) Non observable model; b) Final identified IPN model

Fig 6 .

 6 Fig 6. Two different components running concurrently

  Definition 8. Two transitions t a , t b are called transitions in a two-cycle if S contains the subsequence t a t b t a or the subsequence t b t a t b . The two-cycle transitions set TC of S is given by TC={(t a ,t b )|t a , t b are in a two-cycle}. Example 2. In the sequence S = t 1 t 2 t 3 t 4 t 1 t 2 t 4 t 3 t 5 t 6 t 7 t 4 t 1 t 2 t 3 t 4 t 5 t 6 t 7 t 4 t 1 t 2 t 3 t 4 t 1 t 2 t 3 t 4 t 5 t 6 t 7 t 4 t 5 t 6 t 7 t 4 t 1 t 2 t 3 t 4 t 5 t 6 t 7 t 4 t 5 t 6 t 7 t 4 t 5 t 6 t 7 t 4 t 5 t 6 t 7 t 4 t 5 t 6 t 7 t 4 t 1 t 2 t 4 t 3 t 1 t 2 t 3 t 4 t 1 t 2 t 3 t 4 t 1 t 2 t 4 t 3 t 1 t 2 t 3 t 4 t 1 from Example 1, it can be easily determined Seq={(t 1 ,t 2 ), (t 2 ,t 3 ), (t 3 ,t 4 ), (t 4 ,t 1 ), (t 4 ,t 5 ), (t 5 ,t 6 ), (t 6 ,t 7 ), (t 7 ,t 4 ), (t 2 ,t 4 ), (t 4 ,t 3 ), (t 3 ,t 1 )}. Also, the following precedence relationships are found t 2 ∠t 1 , t 3 ∠t 1 , t 4 ∠t 1 , thus SP(t 1 )={t 1 , t 2 , t 3 , t 4 }. Notice that SP(t j ) is the set of transitions that must invariantly occur to fire t j repeatedly. The rest of the SP sets are SP(t 2 )={t 1 , t 2 , t 3 , t 4 }, SP(t 3 ) = {t 1 , t 2 , t 3 }, SP(t 4 ) = {t 4 }, SP(t 5 ) = {t 4 , t 5 , t 6 , t 7 }, SP(t 6 ) = {t 4 , t 5 , t 6 , t 7 }, SP(t 7 ) = {t 4 , t 5 , t 6 , t 7 }. The set of transitions in a two-cycle is TC= ∅.

  If (t a ,t b ) ∈TC, then there must exist in N a place from t a to t b and a place from t b to t a .Notice that when two transitions are observed consecutively and one is systematically preceded by the other, a causal relationship is found. Also, when two transitions are involved in a two-cycle relation, they are in a causal relationship each other. Observe that all of these relationships are structural, and thus do not depend of the initial marking of the net.Definition 10. The causal relationship set CausalR keeps track of all the causal relationships between transitions in S. CausalR ={(t a ,t b ) | (t a <t b ) and (t a ∠t b or t b ∠t a or (t a ,t b )∈TC)}. From the SP and Seq sets of Example 2, we compute (t 1 ,t 2 )∈CausalR because t 1 <t 2 and t 1 ∠t 2 (Proposition 2). So, CausalR={(t 1 ,t 2 ), (t 2 ,t 3 ), (t 4 ,t 1 ), (t 4 ,t 5 ), (t 5 ,t 6 ), (t 6 ,t 7 ), (t 7 ,t 4 ), (t 2 ,t 4 ), (t 3 ,t 1 )}.

	Example 3.

Proposition 1. If t a ∠t b (t a ∈SP(t b )) then, there exist in N a simple elementary (SE) circuit to which both t a and t b belong. Proposition 2. If t a < t b and either t a ∠t b or t b ∠t a then there must exist in N a place from t a to t b . Proposition 3.

Proposition 4 .

 4 Let t a , t b be two transitions which have been observed consecutively in a complete sequence S in both orders, i.e. (t a , t b )∈Seq, (t b , t a )∈Seq. Then (t a , t b )∉CausalR and (t b , t a )∉CausalR if and only if t a ||t b .

Definition 11.

  Two transitions t a and t b are Sequentially Independent if (t a ,t b ) ∉ CausalR and (t b ,t a ) ∉ CausalR. By inspection of the set CausalR of Example 3, it can be found that t 1 and t 5 are sequentially independent, as well as t 2 and t 6 , and t 3 and t 4 . Proposition 5. Let t a and t b be two transitions in S which have been observed consecutively in both orders (t a <t b and t b <t a ). t a ||t b if: a) t a and t b are Sequentially Independent b) |SP(t a )| >1 and |SP(t b )| >1

  7. Let be two transitions t a and t b which have been observed consecutively in both orders (t a < t b and t b < t a ). t a ||t b if t a and t b are: a) Sequentially Independent, and b) ∃ t k such that t k ∈SP(t b ), t k ∉SP(t a ), and (t a, t k )∈Seq Let be t a and t b two transitions which have been observed consecutively in both orders (t a < t b and t b < t a ). t a ||t b if t a and t b are: a) Sequentially Independent, and b

	Definition 12. The Inverse Systematic Precedence set of a transition SP -1 :T → 2 T contains the transitions which are
	dependent of a common transition to re-enable their firing:
	SP	1	( t	j	)	=	{ t	k	|	t	k	≠	t	j	and	t	j	∈	( t SP	k	)}

-Proposition 8.

  t a , t b ]. In general a place p can be denoted as [t a1 t a2… t al , t b1 t b2… t bh ], where t ai are the input transitions of p and t bi are the output transitions of p, and l=| • p|, h=|p • |. The same place could be used to relate several consecutive transitions. If a transition t k has been observed followed by two transitions t ai , t aj in S (t k <t ai and t k <t aj ), there are two cases to represent such observations into the PN model: a) the case of selection, where t ai , t aj are neither concurrent nor consecutive; thus they are represented with the same place [t k , t ai t aj ], or b) the case of concurrency, where t ai , t aj are concurrent or consecutive and thus they are represented with different places [t k , t ai ] [t k , t aj ]. The same reasoning can be applied to couples of places sharing the same output transition [t ai , t k ] and [t aj , t k ]. In the current example, from CausalR, the causal dependencies [t 4 , t 1 ] and [t 4 , t 5 ] can be represented by a single place [t 4 , t 1 t 5 ] in a selection structure, and [t 2 , t 3 ] and [t 2 , t 4 ] are represented by two places with a single input transition t 2 .

Proposition 9 .

 9 If the IPN model has been correctly build, every computed non-observable place p in N must fulfill the place input/output flow equation: ) is the number of occurrences of t i in S. If there exists a place p such that | • p|=1, then ∀t j ∈ p • , t k ∈ SP(t j ), where t k is the input transition of p. Also, if there exists a place p such that |p • |=1, then ∀t j ∈ • p, t k ∈ SP(t j ), where t k is the output transition of p. Correction rule. If the input/output flow equation or the conditions in Proposition 9 are not satisfied by some place, the arcs relating transitions which are not in CausalR are removed. If there are not pairs in CausalR representing such a place, it is deleted. For the model of fig. 5.b the flow equation is verified for all the transitions.

	( t Occ p ∑ • t i ∈	i	)	( t i Occ • = ∑ ∈ p t i	)	±	1
	where Occ(t i Proposition 10.

Characteristics of Algorithm 1.In step

  Input: The sequence S Output: Non-observable model representing S 1. Compute Seq, SP and TC from S 2. From Seq, SP and TC compute CausalR 3. From Seq and CausalR, compute ConcR 4. Merge transitions as specified in section IV.D 5. Validate and correct places as specified in section IV.E 4 the pairs (t a , t b ) ∈ CausalR are represented by places [t a , t b ] and then merged to build the PN structure; this assures that many pairs in Seq are taken into account. Furthermore, the pairs (t a , t b ) ∈ ConcR, imply that both (t a , t b ) and (t b , t a ) in Seq are represented in the behavior of the PN. Therefore, if Seq'=((Seq \ CausalR) \ ConcR) \ ConcR -1 = ∅ then S can be reproduced by the model built.If Seq'≠∅, some (t b , t a ) in Seq have not been characterized as concurrent or causal. If they are actually sequential, all the verification rules are satisfied; otherwise, they are concurrent and they are corrected in step 5. Once they are corrected, the remaining places relate sequential transitions and thus the sequence S is reproducible.All the procedures and tests derived from the defined sets and propositions involve simple operations on arrays of size |S| or |T|×|T|, which are performed efficiently.

					p 3
			t	1		t	5
				p 1	p 4
				t	2	t	6
	p 2				p 6		p 5
	t	3	a)	t	4	t	7

  4 t 1 t 2 t 3 t 1 t 5 t 6 t 7 t 2 t 3 t 4 t 5 t 6 t 1 t 2 t 3 t 7 t 1 t 4 t 2 t 5 t 4 t 3 t 5 t 1 t 4 t 5 t 2 t 6 t 3 t 1 t 7 t 6 t 7 t 6 t 2 t 3 t 1 t 2 t 3 t 7 t 1 t 4 t 5 t 2 t 4 t 5 t 3 t 4 t 5 t 4 t 5 t 1 t 6 t 7 t 4 t 2 t 3 t 1 t 5 t 6 t 7 t 6 t 2 t 3 t 7 t 4 t 1 t 5 t 2 t 4 t 5 t 6 t 7 t 3 t 1 t 6 t 2 t 7 t 6 t 7 t 6 t 7 t 6 t 3 t 7 t 4 t 1 t 2 t 3 t 5 t 1 t 4 t 5 t 6 t 7 t 2 t 6 t 3 t 1 t 7 t 4 t 2 t 3 t 5 t 6 t 1 t 2 t 3 t 1 t 7 t 4 t 5 t 6 t 7 t 4 t 2 t 3 t 1 t 5 t 4 t 2 t 5 t 6 t 3 t 7 t 6 t 1 t 7 t 4 t 2 t 3 t 5 t 4 t 5 t 1 t 6 t 2 t 7 t 4 t 5 t 6 t 3 t 7 t 1 t 6 t 7 t 2 t 3 t 1 t 2 t 6 t 7 t 4 t 5 t 4 t 5 t 3 t 4 t 5 t 4 t 1 t 5 t 4 t 5 t 4 t 5 t 2 t 6 t 3 t 1 t 2 t 3 t 1 t 7 t 2 t 3 t 1 t 4 t 5 t 6 t 2 t 3 t 1 t 2 t 3 t 1 t 2 t 3 t 1 t 7 t 2 t 6 t 7 t 3 t 1 t 2 t 4 t 5 t 4 yields a PN model composed of two independent subnets shown in Fig. 6.