Joffray Guillory
email: joffray.guillory@cnam.fr.

Daniel Truong

Jean-Pierre Wallerand

Uncertainty assessment of a prototype of multilateration coordinate measurement system

Keywords: dimensional metrology, Large Volume Metrology, Absolute Distance Metre, coordinate measurement system, multilateration technique with self-calibration

Large Volume Metrology is essential to many high value industries to go towards the factory of the future, but also to many science facilities for fine alignment of large structures. In this context, we have developed a multilateration coordinate measurement system, traceable to SI metre, and suitable for outdoor measurements or industrial environments. It is based on a high accuracy absolute distance meter developed in-house and shared between several measurement heads by fibre-optic links. Thus, from these measurement stations, multiple distance measurements of several positions of a target can be performed. At the end, coordinates of the heads and of the different target locations are determined using a multilateration algorithm with self-calibration. In this paper, the uncertainty of this multilateration coordinate measurement system is determined with a consistent metrological approach. First, 13 different sources of errors are listed and quantified. Then, thanks to Monte Carlo simulations, the standard uncertainty on a single absolute distance measurement is assessed to 4.7 µm. This includes the uncertainty contribution of the telemetric system itself, but also the contributions of the mechanical designs of the measurement heads and the target. Lastly, measurements of three-dimensional coordinates of target positions are performed in a control environment, then in a large workshop without temperature control: these measurements validate the uncertainty assessment of the system.

Introduction

Machining and assembly of large components, with measurements carried out as part of the process, is a request in many high value industries, for instance in the automotive [START_REF] Kiraci | Moving towards in-line metrology: Evaluation of a laser radar system for in-line dimensional inspection for automotive assembly systems[END_REF] and aerospace [START_REF] Martin | Metrology enhanced tooling for aerospace (META): A live fixturing, wing box assembly case study[END_REF] fields. In the latter, the assembly tolerances for aircraft wings are lower than 300 μm [START_REF] Martin | Metrology enhanced tooling for aerospace (META): A live fixturing, wing box assembly case study[END_REF]. In this context, large volume metrology (LVM) is the key to enable more accurate coordinate measurements, and thus improve the machine tools, the industrial robots and inspection systems employed. Besides that, the particle accelerators, which deliver beams that become increasingly smaller and more energetic, require positioning even more accurate, at micrometre levels in a range of hundreds of meters [START_REF] Arpaia | Editorial for the special feature on metrology for particle accelerators[END_REF][START_REF] Leão | Engineering survey planning for the alignment of a particle accelerator: part I. Proposition of an assessment method[END_REF]. Indeed, their performances highly depend on the relative alignment between their components. Lastly, LVM can also be applied over radio telescopes used in geodetic and astrometric very-long-baseline interferometry (VLBI) for monitoring of their spatial position and for estimation of their geometrical reference point [START_REF] Lösler | A Modified Approach for Process-Integrated Reference Point Determination[END_REF][START_REF] Lösler | Gravitational deformation of ring-focus antennas for VGOS: first investigations at the Onsala twin telescopes project[END_REF]. In this case, a positioning accuracy better than 1 mm is anticipated in the coming years.

These few examples show the diversity of applications and the levels of uncertainty required. To address these challenges, several instruments commercially available offer sub-millimetre uncertainties for coordinate measurements. For instance, indoor global positioning systems (iGPS) have a level of uncertainty between 85 µm and 500 µm (coverage factor k=1) according to the system configuration [START_REF] Schmitt | Performance Evaluation of iGPS for Industrial Applications[END_REF][START_REF] Muelaner | Verification of the indoor GPS system, by comparison with calibrated coordinates and by angular reference[END_REF][START_REF] Hughes | iGPS Capability Study[END_REF]. As an alternative, photogrammetric systems are also used: some of them, based on dual camera, can achieve uncertainties as low as 50 μm (k=1) [START_REF]Geodetic Systems, Inca 4 -Picture Perfect Measurements -Product Data Sheet[END_REF][START_REF] Martin | Comparative performance between two photogrammetric systems and a reference laser tracker network for large-volume industrial measurement[END_REF]. And when these uncertainty levels become too high for some operations, laser trackers may be adopted. Thanks to a laser distance measurement with two angle measurements, they can provide the threedimensional coordinate of a retroreflector. Manufacturers claim volumetric accuracy of 10 μm + 5 μm/m [12], a value expressed as maximum permissible error (MPE) and which corresponds to an interferometric distance measurement.

In practice, all these instruments present limitations for measurements in industrial environments such as a factory workshop without a temperature control: they are sensitive to temperature gradients. The effect is negligible in terms of distance travelled by the laser beam, but not in terms of beam deflection, up to several tens of micrometres as shown in [START_REF] Pérez Muñoz | Analysis of the initial thermal stabilization and air turbulences effects on Laser Tracker measurements[END_REF][START_REF] Robson | Optimized multi-camera systems for dimensional control in factory environments[END_REF], which has a direct impact on the positioning measurements. Outdoors, this is the air turbulences that could lead to losses of the beam signal, and so to the inability to perform interferometric measurements with laser trackers.

In such cases, the multilateration coordinate measurement technique appears as a good solution since it requires only distance observations of several positions from multiple stations, generally four. This technique has been demonstrated using fringe counting interferometry with commercial laser tracer devices [15] or using frequency scanning interferometry (FSI) systems [START_REF] Dale | Multi-channel absolute distance measurement system with sub ppm-accuracy and 20 m range using frequency scanning interferometry and gas absorption cells[END_REF].

In this paper, we estimate the uncertainty of a multilateration system based on absolute distance measurements, which provides a better flexibility of use and a greater robustness than interferometric systems. Its design lies in the use of a single non-interferometric telemetric system shared between four measurement heads thanks to a network of optical fibres. The adoption of standard optoelectronic and fibre-optic components coming from the telecommunication world and largely available reduces the cost of the system.

The developed non-interferometric telemetric system is based on the measurement of the phase accumulated by a Radio Frequency (RF) modulated light during its propagation in air. Compared to interferometric system, it has a higher synthetic wavelength (tens of millimeters versus hundreds of nanometers) easily traceable thanks to a microwave reference such as an atomic Rubidium clock. Thereby, the system is resilient to wave-front distortions: this reduces the constraints on the optical components (i.e. their surfaces), and provides robustness to atmospheric turbulences due to air stream. Moreover, the system works over a wide range of received power, which allows a range of operation equal to at least 140 m. This is possible thanks to an original design of motorized two-stage rotation system which enables both a full rotation of 360° and a high resolution able to aim at a target positioned 140 m away.

A system based on the measurement of the phase accumulated by a RF modulated light is intrinsically less accurate than interferometric systems, but the micrometric level uncertainty it can reach is sufficient for a wide range of applications. In industrial environment, the limiting factors are the vibrations of the machines, which can involve movements up to 4 µm at maximum velocities of 1 mm/s [START_REF] Muelaner | Absolute multilateration between spheres[END_REF], or the knowledge of the air refractive index, which can be affected by vertical thermal gradient of 1°C/m [START_REF] Muelaner | Large Volume Metrology Technologies for the Light Controlled Factory[END_REF].

Thanks to a micrometric accuracy determined with a consistent metrological approach and traceable to SI metre (International System of units), this system could be used as a reference for on-site calibration of less accurate devices (indoor Global Positioning System, photogrammetric system) or used for determination of the geometrical reference point of radio telescopes.

In this paper, section II presents the multilateration coordinate measurement system that has been developed, then section III lists its different sources of error. Each of them are quantified so that, at the end, the uncertainty on a distance measurement is assessed thanks to Monte Carlo simulations. Lastly, in section IV, the multilateration system is validated experimentally by measurements of three-dimensional coordinates of target positions.

Principle and implementation of the multilateration system

The objective is to determine the three-dimensional positions of targets using absolute distance meters and the multilateration coordinate measurement technique.

As depicted in Fig. 1, a unique telemetric system feeds four measurement heads thanks to a network of optical fibers. By this way, four optical distances can be measured between each head and a target, i.e. a retroreflector. By repeating these measurements for several target positions, an equation system with more observations than unknowns is created. It is therefore possible to determine the coordinates of each target position, but also the coordinates of the four measurement heads: this is the multilateration technique with self-calibration. The telemetric system we have developed is the key element of the system that performs the absolute distance measurements. A laser beam at 1550 nm, propagated through a singlemode fiber, is amplitude modulated by a RF carrier fRF at 4895 MHz, amplified, and sent out towards a 1x4 optical switch which selects one of the four measurement heads, and so the path by which the target will be reached. At the level of the measurement head, the laser beam is emitted in free space and collimated by an off-axis parabolic mirror. The beam diameter at 1% power level is equal to 9.3 mm. After propagation up to the target, a hollow corner cube, the modulated beam returns by the same path to the telemetric system where it is photodetected and where a phasemeter, designed in-house and based on the open-source hardware Red Pitaya, measures the phase shift Φ between this signal and a reference signal directly supplied by the RF synthesizer. This phase shift, measured at an intermediate frequency of 10.75 MHz thanks to a down-conversion by a local oscillator, is directionally proportional to the distance L to measure:

L = (1/2) × ( Φ/(2π) + k ) × ( c/(n×fRF)) (1)
with c the speed of light, n the air group refractive index and k a positive integer. The length measurement is based on the knowledge of the synthetic wavelength, c/f = 61.2 mm, with traceability back to time standards ensured by a reference rubidium (Rb) clock at 10 MHz.

As just explained, the measurement heads are connected to the telemetric system by an optical switch. The four distance measurements are therefore performed sequentially, i.e. one after the other. The developed multilateration system is only working with static targets.

To aim at the target, motorized measurement heads have been developed. They consist of a main gimbal mechanism driven by stepper motors. The platform from where the laser beam emerges can thus rotate around two rotation axes, called standing and transit, with an angular resolution of 440 μrad. However, because such a resolution is insufficient for distances of a few meters, a second gimbal mechanism based on piezoelectric actuators has been implemented. As shown in Fig. 2, it is set up into the main gimbal mechanism: the 90° angle of the offaxis parabolic mirror used for the collimation of the laser beam at the fiber end can be finely tuned by two piezoelectric actuators. Lastly, the target, a hollow corner cube with an aperture of 30°, is also a gimbal mechanism driven by stepper motors that turns towards the head selected by the optical switch for each distance measurement. It has to be noted that our system is compatible with Spherically Mounted Retroreflectors (SMR) usually used by laser trackers. In this case, the SMR has to be turned towards the selected head manually by an operator.

Assessment of the uncertainty on a distance measurement

The assessment of the uncertainty on a distance measurement is quite complex as it depends on many parameters. This task has therefore been realized is two main steps: first, the different sources of errors due to the telemetric system have been identified [START_REF] Guillory | Uncertainty Assessment of Optical Distance Measurements Micrometer Level Accuracy for Long-Range Applications[END_REF], then, the different sources of errors due to the mechanical designs of the measurement heads and the target have been studied [START_REF] Guillory | Assessment of the mechanical errors of a prototype of an optical multilateration system[END_REF]. In this section, the results published in the previous articles [START_REF] Guillory | Uncertainty Assessment of Optical Distance Measurements Micrometer Level Accuracy for Long-Range Applications[END_REF][START_REF] Guillory | Assessment of the mechanical errors of a prototype of an optical multilateration system[END_REF] are summarized, corrected as a result of changes of components in the telemetric system, complemented by new measurements for the mechanical designs of the measurement heads, and merged together to obtain the final value of the uncertainty on a distance measurement. Moreover, the determination of the air group refractive index is also considered.

Sources of errors due to the telemetric system

As shown in formula [START_REF] Kiraci | Moving towards in-line metrology: Evaluation of a laser radar system for in-line dimensional inspection for automotive assembly systems[END_REF], the calculation of the distance is based on the measurement of a phase shift Φ and on the knowledge of the modulation frequency fRF.

The modulation frequency is generated by a RF synthesizer locked on a 10 MHz Rb atomic clock from Microsemi, model SA.31m. As an error in its value leads to a scale error in the distance measurement, its frequency is annually calibrated thanks to a GPS disciplined Rb atomic clock, and at the end, it is only limited by its relative aging rate of ± 1.5 10 -9 per year. Thus, the uncertainty contribution of the modulation frequency is equal to 1.5 10 -9 × L (k=1), i.e. 30 nm for distances of 20 m.

Then, the uncertainty on the phase measurement depends on three sources of errors: the crosstalk, the amplitude to phase coupling, and the random noise.

The crosstalk refers to the addition of a spurious signal at fRF to the ideal measurement signal. To reduce this effect in the developed system, the RF leakages from the emission stages to the reception ones have been minimized thanks to the careful implementation of optical and electromagnetic isolations in some key components. Thus, the Signal to Crosstalk Ratio (SCR) is better than 65 dB, which induces a cyclic error on the distance measurements of amplitude lower than 2.7 µm. With an arcsine distribution, this leads to an uncertainty contribution of the crosstalk of 1.9 µm (k=1).

The amplitude to phase coupling refers to the conversion of an intensity variation of the modulated optical signal into a phase variation of the electrical signal generated by the photodetector. This effect, studied in previous works [START_REF] Guillory | Characterization and reduction of the amplitude-to-phase conversion effects in telemetry[END_REF], induces in our system a linear variation of the measured distance of -0.15 µm/dB. Assuming a worst case indoors, for instance an environment with peak-to-peak power variations up to 5 dB characterized by a uniform distribution, the uncertainty contribution of the amplitude to phase coupling is equal to 0.2 µm (k=1).

Lastly, the random noise refers to the noise performance of the phasemeter with input signals that have been propagated through the whole system. The standard deviation on a phase measurement, in a quiet environment without amplitude variation (no amplitude to phase coupling) and without distance variation (no crosstalk), is 0.17 mrad for 10 ms integration time. Converted into distance, this corresponds to an uncertainty contribution of 0.8 µm (k=1).

Table I I, the dominant source of error is the crosstalk. The criteria to validate a distance measurement is therefore a SCR higher than 65 dB. We have tested the range of operation using one of the measurement head, indoors, in a long corridor. For distances of 70 m and 138 m, the power of the received RF signal used for the phase measurement was around -5 dBm and the crosstalk level vary between -80 dBm and -70 dBm, i.e. a SCR better than 65 dB. The range of operation of our ADM is therefore at least equal to 138 m. Fig. 3 depicts these two distance measurements.

Sources of errors due to the determination of the air refractive index

To deduce a geometric distance from an optical path, the air group refractive index through which the optical beam is propagated should be properly determined. It is generally calculated from an update of the Edlén's formula [START_REF] Bönsch | Measurement of the refractive index of air and comparison with modified Edlén's formulae[END_REF], which depends on the air temperature, the atmospheric pressure, the partial pressure of water vapor, and the CO2 content. In practice, these environmental parameters are measured using local sensors: Table II shows the systematic deviations on a geometric distance due to errors in the measured environmental parameters.

Table II. Sources of errors in the air group refractive index determination based on Edlén's formula. The temperature of air viewed by the optical beam is the most difficult parameter to measure, but also the most critical measurand in the air refractive index determination. Indeed, it should be determined with an accuracy of 1 °C in order to know the geometric distance with an accuracy of about 1 µm/m. For instance, the distance variations observed in Fig. 3, up to 10 µm for a distance of 70 m and up to 18 µm for a distance of 138 m, can be explained by fluctuations of the average temperature along the optical path of 0.15 °C.

Parameters

Sources of errors due to the mechanical designs

Like any mechanical system, the gimbal mechanisms of the measurement heads and the target are not perfectly machined and assembled. This induces misalignments, i.e. additional sources of errors on the geometric distances. These errors are modelled thanks to three parameters depicted in Fig. 4: the beam offset, the beam tilt, and the transit offset. These errors, shared by other instruments such as the laser trackers, total stations and theodolites, are described in literature [START_REF] Muralikrishnan | ASME B89.4.19 Performance Evaluation Tests and Geometric Misalignments in Laser Trackers[END_REF]. There are other mechanical sources of errors in these systems, such as the transit axis tilt (lack of orthogonality between the transit axis and the standing axis), but in this study, we have only included the errors that affect the measured distances (not the angles), which are the only relevant errors of a multilateration system.

The beam offset refers to a laser beam that does not pass by the intersection of the two rotation axes. It is defined by the distances from the rotation axes: a and c for the measurement head, a' and c' for the target. For the target, an offset b' is also considered as in practice the corner cube is not necessarily well aligned with the optical axis. The beam tilt refers to a laser beam that is tilted from its ideal path. It only applies to the measurement head and it Std = 3.6 µm Std = 2.8 µm is defined by two angles: α and γ. Lastly, the transit offset refers to rotation axes that do not cross each other. It is defined by the distance between the transit axis and the standing axis: Toh for the measurement head, Tot for the target. The procedures to measure these nine parameters (a, c, a', b', c', α, γ, Toh and Tot) have been described in [START_REF] Guillory | Assessment of the mechanical errors of a prototype of an optical multilateration system[END_REF] for the target and for an earlier design of the measurement heads not including a piezoelectric gimbal mount. The addition of a piezoelectric gimbal mount has not changed the way to measure these parameters, which is based on a double-centering method (also called double-face measurement) as detailed in [START_REF] Guillory | Assessment of the mechanical errors of a prototype of an optical multilateration system[END_REF]. This method consists in pointing the same object under two different orientations of the head (or of the corner cube), thanks to a first rotation of 180° around the standing axis and then a second 180° rotation around the transit axis. If the optical beam does not emerge from (or if the corner cube is not positioned at) the intersection of the two rotating axes, orthogonally to the plane formed by these axes, this procedure will produce two measurements with opposite errors from which the beam offset, the beam tilt, and the transit offset could be extracted.

With the piezoelectric gimbal, the beam offset and the beam tilt depend on the values of the piezoelectric actuators. This is easily understandable for the beam tilt because the variations of the angles α and γ correspond actually to the fine adjustment of the viewing angles of the measurement heads. For each measurement head, the beam offset and the beam tilt parameters have been measured for the following typical cases: when the two piezoelectric actuators are in their zero-position (0 V), and when they are driven in half (75 V) and in full travel-range (150 V). Results are depicted in Fig. 5 for the beam offsets and in Fig. 6 for the beam tilts. As depicted in Fig. 5, the beam offset can vary widely depending on head, up to 285 µm for the measurement head C. Such a value cannot be explained by the design of the piezoelectric gimbal mount depicted in Fig. 2, and the 20 µm stroke of its piezoelectric actuators. The explanation probably comes from the way the position of the laser beam is experimentally determined: it is defined as the position where the maximum amount of optical power is detected. The latter may change, not because the whole laser beam moves, but because the power distribution within the laser beam evolves with the tilt applies to the parabolic mirror. In other words, the distortions of the gaussian beam due to variations of the angle of incidence on the parabolic mirror are interpreted as changes of the beam offset.

The different measured values of the beam offset have been plotted in Fig. 7 as a histogram, taking into account all the measurements made on the four heads and for different voltages applied to the piezo actuators. We recognize a uniform distribution. Therefore, for simplicity reasons, the beam offset of all the measurement heads is finally consider as a random variable characterized by a unique uniform distribution on an interval between -320 µm and + 320 µm. The value of 320 µm corresponds to the worst beam offset observed in Concerning the beam tilt, which corresponds to the fine adjustment of the viewing angles of the measurement heads, it follows a uniform distribution for each head since all the positions of the piezoelectric actuators can be selected with equal probability. However, the range slightly depends on the measurement head and on the beam tilt axis. Uncertainties on the beam tilts measurements (i.e. uncertainty bars in Fig. 6) have been estimated, in worst case, to 0.01°.

At the end, Table III summarizes the results for each head, with the values measured for each misalignment parameter (or the range of possible values), and their probability distributions. In Tables III and IV, the contributions of the error parameters to a distance measurement have not been stated. Nevertheless, they are detailed in [START_REF] Guillory | Assessment of the mechanical errors of a prototype of an optical multilateration system[END_REF]. In short, their contribution can be low when considering separately these different sources of error. For instance, the beam offsets a and c (or a' and c') are error components perpendicular to the measured distance, which not produce significant errors. Concerning the beam tilts α and γ, when they are considered independently, they do not induce error on the distance measurement. However, the transit offset can impact strongly the measurements. Its contribution, which depends on the elevation angle φ, is as follows [START_REF] Guillory | Assessment of the mechanical errors of a prototype of an optical multilateration system[END_REF]: Rm² = Rt² + Toh² + 2 Toh Rt cos(φ) (or Rm² = Rt² + Tot² + 2 Tot Rt cos(φ') for the target) [START_REF] Martin | Metrology enhanced tooling for aerospace (META): A live fixturing, wing box assembly case study[END_REF] with Rm the measured distance and Rt the true distance (R stands for range). This is the main error contribution to a distance measurement, but in practice, it can be corrected thanks to the measurement of the elevation angle.

Assessment of the uncertainty on a distance measurement

In this part, the uncertainty on a distance measurement is quantified by Monte Carlo simulations using the uncertainty distributions of the different sources of errors listed above. However, the contribution of the air group refractive index is not developed as it does not depend on the developed system, but on the accuracy of the local sensors we used and on the fluctuations of the environmental parameters, especially the temperature.

First, the contribution of the telemetric system is quantified. To determine the impact of the four additive sources of errors (Table I) on distance measurements up to 20 m, a Monte Carlo simulation has been performed since each error has a different distribution. Thus, for each source of errors, 200 000 random samples have been generated following the probability distributions defined previously. For a given distance, the measured distance (Rm for measured range) is obtained from the true one (Rt) by:

Rm = Rt + 1.5•10 -9 G Rt + 2.7•10 -6 sin(2π U) + 2 sqrt(3) 0.2•10 -6 U + 0.8•10 -6 G ( 3
)
where G is a random value following a gaussian distribution centered at zero and of variance 1, and U is a random value between -0.5 and +0.5 following a uniform distribution. Thus, the four additive errors have been considered: the length-dependent gaussian error of 1.5 nm/m due to the modulation frequency, the cyclic error of 2.7 µm due to the crosstalk, the uniform error of 0.2 µm due to the amplitude-to-phase coupling, and the gaussian error of 0.8 µm due to the random noise.

After running the Monte Carlo simulation for different distances to measure, from 20 cm to 20.2 m by step of 80 cm, and after the addition of all these errors (26 tested distances × 200 000 random samples = 5.2 millions of values), we obtain the distribution shown in Fig. 9, on the left. It has to be noted that the Monte Carlo simulations made for different distances have not shown a length dependence of the errors. That is why they have been joined in a same distribution. Moreover, in a predictable way, the shape of the distribution demonstrated that the dominant source of error is the crosstalk with its arcsine distribution.

Then, the contribution of the mechanical designs has been studied. In [START_REF] Guillory | Assessment of the mechanical errors of a prototype of an optical multilateration system[END_REF], the true distance has been modelled as a function of the measured one. For the measurement heads, the model takes into account the misalignment parameters of Table III, but also its elevation angle φ. More specifically, the true distance (Rt) is obtained from an intermediate result (Ri) corresponding to a distance measurement when only beam offset and beam tilt are considered, and then, the measured value (Rm) is obtained from this intermediate result (Ri) when the transit offset is also added. The relationship between the true distance and the measured one can so be determined from this system of two equations: Rt² = Ri² + a² + c² -2 a sqrt( Ri² + c² -2 c Ri sin(γ) ) sin(α) -2 c Ri sin(γ) Rm² = Toh² + Ri² + 2 Toh Ri cos(α) cos(φ-γ) (4) For the target, the model considers the misalignment parameters of Table IV, the elevation angle of the retroreflector φ', and also its misalignment with respect to the optical axis θ'' et φ''. The model is even more complex:

Rm² = Rt² + Tot² + a'² + b'² + c'² + 2 a' Rt sin(θ") + 2 b' [ Tot cos(φ') -Rt cos(θ") cos(φ") ]
+ 2 c' [ Rt sin(φ") -Tot sin(φ') ] -2 Tot Rt [ cos(φ') cos(θ") cos(φ") + sin(φ') sin(φ") ]

These two models involve many variables with different distributions. It is consequently difficult to calculate the propagation of the uncertainties to assess their impact on a distance measurement, and again, a Monte Carlo analysis has been preferred. In practice, the contribution of the measurement heads and the one of the target have been studied separately: for each case a Monte Carlo analysis has been performed following the diagram depicted in Fig. 8. First, different distances to measure have been considered, from 20 cm up to 20 m. Then, errors have been applied on these distances following the model we have defined. To this end, 200 000 random samples of the beam offset, beam tilt, and transit offset have been generated following the probability distributions defined previously. In addition to this, other parameters have been considered: for instance, for the target, elevation angle φ' and misalignment angles θ'' et φ'' have to be included in the model. As all angle positions can be achieved with equal probability, a random sampling characterized by a uniform distribution has been used.

In theory, all systematic errors can be corrected, and in that case, we consider the uncertainties of the corrections themselves [START_REF] Mesures | Evaluation of measurement data -An introduction to the "Guide to the expression of uncertainty in measurement" and related documents[END_REF]. In our case, the systematic errors due to beam offsets and beam tilts are not corrected, therefore these bias elements are transformed into variances. errors due to transit offsets are corrected using the defined models, the estimations of their values (the ones in Tables III or IV), and the values of the elevation angles (φ for the measurement head or φ' for the target).

At the end, a set of measured distances affected by mechanical errors are obtained. And because no length dependence of these errors is visible, all the errors have been joined in a same distribution. Thus, the distribution of the errors on the measured distances due to the mechanical design of the measurement head A and of the target is depicted in Fig. 9. The results for the three other measurement heads are similar to the measurement head A. Figure 9. Distribution of the errors on a distance measurement after Monte Carlo analysis. 68% of the errors are within the red part, and 95% within the red and blue parts.

The mechanical errors due to different angular orientations of the rotating reflector, characterized by a gaussian distribution with the standard deviation of 3.9 µm, dominate the other sources of errors. In contrast, the mechanical errors of the measurement heads are much lower, with a standard deviation of only 1.4 µm and a distribution that can be characterized by a Lorentzian.

At the end, assuming no error is brought by the determination of air refractive index, and considering additive errors for the three contributions, i.e. the telemetric system, the mechanical designs of the measurement heads and of the target, we obtain an uncertainty on a distance measurement of 4.7 µm (k=1). Its distribution, depicted in Fig. 10 and obtained by summing the three distributions in Fig. 9, is very close to a gaussian curve due to errors dominated by the mechanical design of the target.

Figure 10. Error distribution on a distance measurement performed by the multilateration system.

Experimental validation

The validation of the developed multilateration system has been done through two experiments, one over a small volume of one cubic meter, and one over a large volume close to an industrial environment. In both cases, the four measurement heads have been set up to form, as far as possible, a regular tetrahedron, the optimal configuration for multilateration algorithm, as explained in [START_REF] Takatsuji | Restriction on the arrangement of laser trackers in laser trilateration[END_REF].

Small volume configuration

In the first case, the four measurement heads have been positioned on an optical table at the coordinates: In this volume of about one cubic meter, 14 target positions have been distributed. Among them, the points 1, 2 and 3 are a triplet of aligned target positions mounted on a same breadboard as depicted Fig. 11. During this experiment, the breadboard has been displaced to also compose the points 6, 7 and 8, then the points 12, 13 and 14. If the developed coordinate multilateration measurement system works properly, we should obtain the same distances between the different positions of this triplet at different locations. Points numbered 4, 5, 9, 10 and 11 are isolated points as depicted in black in Fig. 11.

A = [
In this small volume, the maximal measured distance is 1.47 m. Hence, the piezoelectric gimbal mechanism has not been used in this experiment: the two piezoelectric actuators were remained in their zero-position. In practice, each position is defined by a tribrach, i.e. an adapter able to receive the target. Thus, an operator can easily move the target from one position to another. For each position, four optical distance measurements have been performed, one per head, to obtain at the end a set of 56 measurements. They have been corrected by the air refractive index. In this experiment, in an air-conditioned room, only one environmental sensor has been used: temperature was around 22.5 °C, pressure around 991.0 hPa, and humidity around 36%.

Large volume configuration

The second configuration is a large volume of 5. As previously, 14 target positions have been distributed, with some of them largely outside the volume bounded by the measurement heads. In this way, a wide variety of distances can be measured, from 0.34 m to 11.54 m. Besides, in this volume, there are three couples of target positions mounted on a same breadboard: again, if the developed coordinate multilateration measurement system works properly, we should obtain the same interpoint distance for these three couples. Additionally, as depicted in green in Fig. 12, there are three aligned pillars.

As previously, a set of 56 distances has been measured, which has taken about 3.5 hours. During this time, the temperature, around 19 °C, has increased by 0.6 °C, and the pressure, around 1016.5 hPa, has decreased by 0.7 hPa. Additionally, the sensor at a height of 2.4 m has recorded temperatures 1 °C higher than the other sensors at a height of 0.3 m, 0.7 m and 1.4 m, which indicates the presence of a vertical temperature gradient. For determination of the air refractive index, this gradient has been ignored as most of the targets have a height less than 1.5 m. Nevertheless, we assume an error up to 0.5 °C on the air temperature, i.e. an additional distance error of 5.4 µm for the longer distance.

Data processing

To determine the 14 positions Tj of the target and the 4 positions Hi of the measurement heads using length measurements only, we have developed a multilateration algorithm with self-calibration based on [START_REF] Zhuang | Self-calibration and mirror center offset elimination of a multibeam laser tracking system[END_REF]. It consists in a nonlinear optimization problem where the quadratic sum of the differences between the squared distances measured by our absolute distance meter (di,j²) and the squared distances calculated from the positions provided by algorithm is minimized:

cost function = Σ ( di,j² -|| Hi-Tj ||² )² (6)
The unknown variables to determine are, of course, the coordinates of the 14 target positions and of the 4 measurement heads, but also the instrument offset of each measurement head. The latter is an additive constant that compensates from delays in electrical cables and optical paths. Thanks to these corrections, absolute distance measurements can be achieved between each measurement head and the target.

The convergence of the algorithm depends on the chosen initial values, and as explained in [START_REF] Zhuang | Self-calibration and mirror center offset elimination of a multibeam laser tracking system[END_REF], they have to be sufficiently accurate, in the order of magnitude of 10 mm. In our case, they have been determined by measurements of distances and angles of one of the four measurement heads, in the same way as a laser tracker. The 14 initial positions of the target j are given by: xj = di,j × sin(φi) × cos(θi), yj = di,j × sin(φi) × sin(θi), zj = di,j × cos(φi) [START_REF] Schmitt | Performance Evaluation of iGPS for Industrial Applications[END_REF] with di,j the measured distances, θi and φi the azimuth and elevation angles of the measurement head i. Once the 14 intial target positions have been determined, the initial coordinates of the four measurement heads are calculated by a classical multilateration algorithm [START_REF] Norrdine | An Algebraic Solution to the Multilateration Problem[END_REF].

To minimize the cost function described in formula ( 6) correctly, at least nine target positions have to be measured [START_REF] Zhuang | Self-calibration and mirror center offset elimination of a multibeam laser tracking system[END_REF][START_REF] Chen | Optimization method for solution model of laser tracker multilateration measurement[END_REF]. This is the minimum of points required to have an equation system with more observations than unknowns. In our experiments, with 14 points, this condition is satisfied: this redundancy helps the algorithm to converge.

Experimental results

Fig. 13 presents the results for the two experiments after applying the multilateration algorithm with selfcalibration. In this figure, the error is the difference between the distances measured by the telemetric system and the distances deduced from the multilateration algorithm that provides the coordinates of the measurement heads and of the targets. For each target position, four errors are thus calculated, one per head. At the end, we obtain small differences between the distances measured by our absolute distance meter, and the ones deduced from the positions provided by the multilateration algorithm with self-calibration. Standard deviation on the error is better than 4 µm in both cases, and surprisingly, better results have been obtained over the large volume where the environmental conditions were less stable, possibly because the position of the measurement heads were more stable thanks to the use of real pillars.

Experimental results are consistent with the uncertainty obtained by Monte Carlo simulation on the distances measured by our absolute distance meter, i.e. 4.7 µm (k=1). Nevertheless, in order to ascertain that the coordinates provided by the multilateration algorithm with self-calibration are correct, some distances between different target positions have been verified.

In the small volume of one cubic meter, the relative distances between the different positions of the triplet have been measured in a direct way by the absolute distance meter (it means without multilateration). To this end, the measurement head A has been positioned in the alignment of the three positions of the triplet, then by moving the target from one position to another, we can determine the relative distances between them: The uncertainties on the relative distances correspond to the standard deviation over 9 successive measurements. These uncertainties include performances of the instrument, the misalignment of the target (even if we try to optimize it) and the centering repeatability of each tribrach, i.e. when the target is removed and mounted again in the tribrachs with unchanged levelling and orientation of the carrier. According to [START_REF] Braun | Suppression of systematic errors of electronic distance meters for measurement of short distances[END_REF], centering repeatability is between 2 and 5 µm.

In the large volume, similarly, the interpoint distances of the three couples of target positions, equal to 324.071 mm, and of the three aligned pillars, equal to 2 232.037 mm, 6 700.471 mm and 8 932.509 mm, have also been measured in a direct way. Comparisons of these direct distance measurements (reference values with uncertainties around 4 µm) with the distances calculated from coordinates obtained by the multilateration algorithm are presented in Tables V and VI. In small volume configuration (Table V), the distance differences for the different positions of the triplet are below 10 µm in 8 out of 9 cases, with standard deviation on these differences of 8.0 µm. This is compatible with the combined standard uncertainty between the telemetric system itself (4.7 µm at k=1) and the comparison itself (4 µm at k=1 for the reference distances).

For large volume configuration (Table VI), one of the couple comparisons leads to 24 µm of error. Otherwise all the comparisons between multilateration measurements and direct measurements are within the uncertainty bars of coverage factor k=2 (95% confidence). Nevertheless, the large volume configuration leads to a standard deviation on the differences a little bit higher, 12.8 µm, possibly due to bad estimates of the air refractive index.

In these comparisons between direct distance measurements and distances calculated from coordinates obtained by multilateration, the uncertainty induced by the multilateration algorithm with self-calibration has been assumed to be zero. In practice, and as demonstrated by the experimental results, the multilateration algorithm with selfcalibration does not induce a significant additional uncertainty compared to the uncertainty induced by the developed system itself.

The uncertainty induced by the multilateration algorithm is well documented in literature for classic multilateration problems through the concept of dilution of precision (DOP) used in global-navigation-satellite-system (GNSS) positioning [START_REF] Langley | Dilution of precision[END_REF]. In the case of the multilateration with self-calibration, more complex, the uncertainty is difficult to assess by a mathematical approach and simulations are preferably used. Thus, we have generated 4 positions of measurement heads and 14 positions of targets with respect to the previous setups depicted Fig. 11 and12. From these points, we have calculated the theoretical distances between the measurement heads and the targets, then added to them the instrument offsets and a random noise following a gaussian distribution centered at zero and of standard deviation 4.7 µm. After running our multilateration algorithm with self-calibration with these noisy values, we have obtained errors on the interpoint distances similar to the ones shown in Tables V and VI. Over several simulations, the standard deviations on these errors can vary from 2.0 µm to 10.7 µm (average value of 6.6 µm over 50 simulations) for the small volume configuration, and from 2.2 µm to 9.2 µm (average value of 5.4 µm) for the large volume configuration, which is compatible with our results. This aspect, i.e. the uncertainty induced by the multilateration algorithm, would deserve a more extensive study in a detailed paper.

Conclusion

A reference multilateration system with metrological traceability to SI metre has been realized. It is composed of a telemetric system, four measurement heads, and a target and it has been studied with a consistent metrological approach. The uncertainty contribution of the telemetric system itself, an absolute distance meter based on the phase shift measurement of an intensity-modulated light, is equal to 2.1 µm (k=1). The uncertainty contribution of misalignments in the gimbal mechanisms of each measurement head and the target is 1.4 µm and 3.9 µm, respectively (k=1). This leads at the end to an uncertainty on a distance measurement equal to 4.7 µm (k=1).

The developed multilateration system has been then validated through two experiments, one over a small volume of one cubic meter, and one over a large volume close to an industrial environment. In both cases, the threedimensional coordinates of 14 target positions have been determined: the multilateration algorithm with selfcalibration has perfectly converged and standard deviations on the errors lower than 4 µm have been demonstrated.

In parallel, distances between some of these target positions have been measured in a direct way by the absolute distance meter. The comparison between these direct distance measurements, with uncertainties around 4 µm, and the distances calculated from coordinates obtained by the multilateration system shows a good consistency. The residuals, for interpoint distances up to 8.9 m, showed that the standard deviation on the differences between the two systems was better than 13 µm with a maximum observed difference of 24 µm.

The developed multilateration coordinate measurement system, which offers high accuracy at a reasonable cost, is therefore a promising solution for large volume measurement, especially in uncontrolled environments or outdoors thanks to non-interferometric absolute distance measurements. However, there are still ongoing developments to improve the proposed system.

In a close future, the target could be replaced by a retroreflecting sphere of glass refractive index n=2. Such a solution is simpler to use than a corner cube mounted a gimbal mechanism thanks to an acceptance angle as large as 360° in principle. The latter has therefore some advantages, for instance it induces no mechanical source of error and it can be seen by all the measurement heads simultaneously. However, according to preliminary measurements, it presents a bad optical reflection with an optical signal reduced by at least 18 dB for short distances, and more for longer distances. This will limit the range of our system: first results have demonstrated a possible range of 20 m with an appropriate diameter of the sphere.

Until now, the alignments between the measurement heads and the target were performed manually thanks to a remote control of the motors of the gimbal mechanisms by Wi-Fi. However, we are working on a procedure for automated alignment and measurement.

Lastly, we plan to perform a direct comparison of the three-dimensional coordinates obtained by our multilateration system with the ones obtained by an interferometric laser tracker. Thus, the system could be validated by an independent calibrated reference measurement system. Figure 9. Distribution of the errors on a distance measurement after Monte Carlo analysis. 68% of the errors are within the red part, and 95% within the red and blue parts. Figure 10. Error distribution on a distance measurement performed by the multilateration system. Figure 11. Layout of the 4 measurement heads and 14 target positions, with a photograph of the optical table on the left, the triplet of aligned target positions on the top, and a graph on the right (top view). Figure 12. Top view of the workshop configuration. Figure 13. Experimental results for the small volume configuration on the left, and for the large volume configuration on the right.

Figure Captions

Table Captions

Table I. Sources of errors of the telemetric system.

Table II. Sources of errors in the air group refractive index determination based on Edlén's formula.

Table III. Sources of errors of the four measurement heads.

Table IV. Sources of errors of the target. Table V. Comparison for the small volume configuration. Table VI. Comparison for the large volume configuration.

Figure. 1 .

 1 Figure. 1. Photograph of the developed multilateration system composed of a shared telemetric system, four measurement heads, and one target.

Figure 2 .

 2 Figure 2. 3D model of the motorized measurement heads and the target, with details on the piezoelectric gimbal mount.

Figure 3 .

 3 Figure 3. Distance measurements over 70 m and 138 m without air group refractive index correction. The measured absolute distances are the sum of the values D in the legend and of the relative distance curves.

Figure 4 .

 4 Figure 4. Illustration of the three misalignment errors in the case of the measurement head.

Figure 5 .

 5 Figure 5. Beam offsets of the heads as a function of the values of the piezoelectric actuators.

Figure 6 .

 6 Figure 6. Beam tilts of the heads as a function of the values of the piezoelectric actuators.

  Fig 5, uncertainty included (k=1).

Figure 7 .

 7 Figure 7. Distribution of the measured beam offsets of the four heads.

  Sources of errors of the four measurement heads. Parameters Description Value for head A Value for head B Value for head C Value for head γ Beam tiltα ∈ [-0.05°,0.12°] γ ∈ [-0.09°,0.02°] α ∈ [0.07°,0.24°] γ ∈ [-0.16°,-0.03°] α ∈ [-0.11°,0.08°] γ ∈ [-0.07°,0.05°]α ∈ [0.09°,0.29°] γ ∈ [-0.08°,0.02°] the misalignment parameters of the target: the beam offset and the transit offset. The latter, measured in[START_REF] Guillory | Assessment of the mechanical errors of a prototype of an optical multilateration system[END_REF], follow a gaussian distribution of standard deviation 2 µm.TableIV. Sources of errors of the target.

Figure 8 .

 8 Figure 8. Diagram of the Monte Carlo analysis.

  504, 938, 491], B = [0, 0, 0], C = [524, -7, 909] and D = [1050, 0, 0] (units in mm), which is a tetrahedron with edges of 1091 mm ± 71 mm.

Figure 11 .

 11 Figure 11. Layout of the 4 measurement heads and 14 target positions, with a photograph of the optical table on the left, the triplet of aligned target positions on the top, and a graph on the right (top view).

  6 m (x) × 10.3 m (y) × 2.6 m (z) where four measurement heads have been positioned at the coordinates: A = [0, 0, 0], B = [2736, 56, 1406], C = [5289, 0, 0] and D = [2692, 6236, 793] (units in mm).

Figure 12 .

 12 Figure 12. Top view of the workshop configuration.

Figure 13 .

 13 Figure 13. Experimental results for the small volume configuration on the left, and for the large volume configuration on the right.

  || P1-P2 ||ref = 150.022 mm ± 4 µm, || P2-P3 ||ref = 174.615 mm ± 4 µm, and || P1-P3 ||ref = 324.636 mm ± 3 µm.

Figure. 1 .

 1 Figure. 1. Photograph of the developed multilateration system composed of a shared telemetric system, four measurement heads, and one target.Figure 2. 3D model of the motorized measurement heads and the target, with details on the piezoelectric gimbal mount. Figure 3. Distance measurements over 70 m and 138 m without air group refractive index correction. The measured absolute distances are the sum of the values D in the legend and of the relative distance curves. Figure. 4. Illustration of the three misalignment errors in the case of the measurement head.Figure 5. Beam offsets of the heads as a function of the values of the piezoelectric actuators.Figure 6. Beam tilts of the heads as a function of the values of the piezoelectric actuators.

Figure 2 .Figure 3 .

 23 Figure. 1. Photograph of the developed multilateration system composed of a shared telemetric system, four measurement heads, and one target.Figure 2. 3D model of the motorized measurement heads and the target, with details on the piezoelectric gimbal mount. Figure 3. Distance measurements over 70 m and 138 m without air group refractive index correction. The measured absolute distances are the sum of the values D in the legend and of the relative distance curves. Figure. 4. Illustration of the three misalignment errors in the case of the measurement head.Figure 5. Beam offsets of the heads as a function of the values of the piezoelectric actuators.Figure 6. Beam tilts of the heads as a function of the values of the piezoelectric actuators.

Figure 5 .

 5 Figure. 1. Photograph of the developed multilateration system composed of a shared telemetric system, four measurement heads, and one target.Figure 2. 3D model of the motorized measurement heads and the target, with details on the piezoelectric gimbal mount. Figure 3. Distance measurements over 70 m and 138 m without air group refractive index correction. The measured absolute distances are the sum of the values D in the legend and of the relative distance curves. Figure. 4. Illustration of the three misalignment errors in the case of the measurement head.Figure 5. Beam offsets of the heads as a function of the values of the piezoelectric actuators.Figure 6. Beam tilts of the heads as a function of the values of the piezoelectric actuators.

Figure 6 .

 6 Figure. 1. Photograph of the developed multilateration system composed of a shared telemetric system, four measurement heads, and one target.Figure 2. 3D model of the motorized measurement heads and the target, with details on the piezoelectric gimbal mount. Figure 3. Distance measurements over 70 m and 138 m without air group refractive index correction. The measured absolute distances are the sum of the values D in the legend and of the relative distance curves. Figure. 4. Illustration of the three misalignment errors in the case of the measurement head.Figure 5. Beam offsets of the heads as a function of the values of the piezoelectric actuators.Figure 6. Beam tilts of the heads as a function of the values of the piezoelectric actuators.

Figure 7 .Figure 8 .

 78 Figure 7. Distribution of the measured beam offsets of the four heads. Figure8. Diagram of the Monte Carlo analysis. Figure9. Distribution of the errors on a distance measurement after Monte Carlo analysis. 68% of the errors are within the red part, and 95% within the red and blue parts. Figure10. Error distribution on a distance measurement performed by the multilateration system. Figure11. Layout of the 4 measurement heads and 14 target positions, with a photograph of the optical table on the left, the triplet of aligned target positions on the top, and a graph on the right (top view). Figure12. Top view of the workshop configuration. Figure13. Experimental results for the small volume configuration on the left, and for the large volume configuration on the right.

  . Sources of errors of the telemetric system.

		Parameter	Description	Value	Contribution	Error distribution
	1	u RF	modulation frequency	fRF = 4895 MHz	1.5 • 10 -9 × L	gaussian distribution
	2 u crosstalk	crosstalk	SCR > 65 dB	1.9 µm	arcsine distribution
	3	u AM/PM	amplitude-to-phase coupling	5 dB power variations	0.2 µm	uniform distribution
	4	u noise	random noise	σφ = 0.17 mrad	0.8 µm	gaussian distribution
	As shown in Table

Table V .

 V Comparison for the small volume configuration.

	Std = 3.5 µm		Std = 1.4 µm
	Volume of 1 m 3	Large Volume
	Triplet || P1-P2 || -|| P1-P2 ||ref || P2-P3 || -|| P2-P3 ||ref || P1-P3 || -|| P1-P3 ||ref
	1	-2.7 µm	2.0 µm	-2.4 µm
	2	-9.8 µm	16.0 µm	4.5 µm
	3	-7.0 µm	8.7 µm	0.0 µm
		Table VI. Comparison for the large volume configuration.
	Couple	|| P1-P3 || -|| P1-P3 ||ref	The three aligned pillars
	1	24.2 µm	|| P1-P2 || -|| P1-P2 ||ref	-9.8 µm
	2	8.6 µm	|| P2-P3 || -|| P2-P3 ||ref	-3.9 µm
	3	-4.9 µm	|| P1-P3 || -|| P1-P3 ||ref	11.3 µm

Acknowledgment

This work was partially funded by Joint Research Projects (JRPs) 17IND03 LaVA and 18SIB01 GeoMetre, projects that have received funding from the European Metrology Programme for Innovation and Research (EMPIR) co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme.


Mariam Bouzid
email: mariam.bouzid@utt.fr

Oussama Masmoudi
email: oussama.masmoudi@utt.fr

Alice Yalaoui
email: alice.yalaoui@utt.fr

Agadir -Morocco "New advances and challenges for sustainable and smart industries" TIME-INDEXED FORMULATIONS FOR THE ORDER ACCEPTANCE SCHEDULING PROBLEM UNDER ENERGY ASPECTS

Keywords: Energy, Order Acceptance Scheduling, Time-indexed formulation

To comply with the new challenges of sustainability, the industrial sector is revising its means of supply and production. This entails, for instance, optimizing energy consumption and costs at the operational level. In this vein, this research presents an Order Acceptance Scheduling problem (OAS) on a single machine under electricity time-of-use tariffs and taxed carbon emission periods. The objective is to maximize the total profit. This problem arises when a company decides to select and process a subset of orders only if it is possible within a predetermined time-window. Therefore, the number of possible schedules grow at a factorial rate. To tackle this NP-hard problem, two time-indexed formulations are provided. Finally, to assess the performance of the proposed models, a comparative analysis against a classical formulation is conducted.

INTRODUCTION

For centuries, industry has been a vector for social and economic prosperity, shimmering an indefinite growth. Although, new concerns about environmental issues, such as resource depletion and greenhouse gases (GHG) emission, cast doubt upon this vision. According to the International Energy Agency, the industrial sector accounted for more than a third of the energy used in the world in 2017, which is responsible for climate changes and jeopardizing social advances at the same time. Industry 4.0through the development of smart decision tools and the modernizing of equipmentallows a rational and an accountable response to the sustainability stakes, which aim at maintaining environmental, social and economic viability. At the strategic level, this is embodied by Corporate Social Responsibility (CSR) that advocate for ethics and green sustainability from global reference frameworks. At the operational level, this includes optimizing energy consumption, costs and carbon footprint. To perform efficient demand management, energy suppliers developed preferential rate, designated as timeof-use (TOU) tariffs, at specific times of the day. In the meantime, industrials must abide by the rules on regulation of GHG emissions, which are reflected by the implementation of carbon emission taxes by governments. We can mention that the objectives for reducing GHG emissions were reinforced at COP22 in Marrakesh, Morocco.

In this context, we present two time-indexed models for the OAS problem with electricity TOU tariffs and taxed carbon emission periods. This research follows up the work of [START_REF] Chen | Order acceptance and scheduling problem with carbon emission reduction and electricity[END_REF] which presented a disjunctive model.

The paper is organised as follows. Section 2 includes a review on OAS problems and scheduling under electricity TOU tariffs. Section 3 states the problem and presents the solution approach. Section 4 features the resolution method. Section 5 presents the computational results. The last section concludes the paper and draws perspectives.

STATE OF THE ART

The standard OAS problem is a double-decision problem that consists in the selection and the sequencing of a subset of ordersamong 𝑛with the objective to maximize the total profit. [START_REF] Slotnick | Order acceptance and scheduling: A taxonomy and review[END_REF] proposes a literature review on this topic, indicating that OAS are studied for both single and multi-machines systems and with various job characteristics such as preemption, release date or setup. These problems are generally known to be NP-hard as demonstrated in [START_REF] Palakiti | Order acceptance and scheduling: overview and complexity results[END_REF]. At worst the number of possible schedules is ∑ 𝑘! 𝑛 𝑘=1 , while in standard scheduling problems, all the orders are accepted and thus only 𝑛! sequenceswhich are all the possible permutations of 𝑛 elements without repetitioncan be obtained. (Oğuz et al., 2010) address the OAS problem with release dates, setup times and time-related penalties using a disjunctive Mixed Integer Linear Program (MILP). For the same problem, [START_REF] Cesaret | A tabu search algorithm for order acceptance and scheduling[END_REF] propose a Tabu Search and [START_REF] Silva | Exact and heuristic algorithms for order acceptance and scheduling with sequence-depend-ent setup times[END_REF] provide an efficient arctime-indexed model. However, in the literature, time-indexed formulations are not as developed for OAS problem.

MOSIM '20 -November 12-14, 2020 -Agadir -Morocco Energy considerations are essential for both economic and environmental reasons. With energy prices increase, demanding specification and taxes, the operational level takes a crucial part in the efforts for sustainability. Pricing policies, especially TOU rate, are largely studied for single and more complex systems. For instance, [START_REF] Aghelinejad | Production scheduling optimisation with machine state and time-dependent energy costs[END_REF] exploit machine states mechanism to minimize total energy costscomprising idle, transition and processing energyon a single machine with a predetermined sequence. The latter use an on/off time-indexed model. [START_REF] Che | An efficient greedy insertion heuristic for energy-conscious single machine scheduling problem under time-of-use electricity tariffs[END_REF] formulate a timeindexed MILP and designed a greedy heuristic for the single machine under TOU electricity tariffs in order to minimize total energy costs. The same authors in [START_REF] Che | Energy-conscious unrelated parallel machine scheduling under timeof-use electricity tariffs[END_REF] investigate the unrelated parallel machine under TOU tariffs with energy costs minimization. [START_REF] Ho | Electricity cost minimisation for optimal makespan solution in flow shop scheduling under time-of-use tariffs[END_REF] jointly enhance energy costs and makespan for a two-machine flow shop under TOU tariffs. Finally, for a job-shop system with TOU tariffs and peak-power considerations, [START_REF] Masmoudi | Job-shop scheduling problem with energy consideration[END_REF] employ a time-indexed formulation for the problem of minimizing energy costs. In the literature, energy aspects are also integrated as constraints. In [START_REF] Liao | Multi-objective optimization of single machine scheduling with energy consumption constraints[END_REF], weighted completion time and weighted tardiness are minimized on a single machine with a periodic threshold on energy consumption. [START_REF] Fang | Flow shop scheduling with peak power consumption constraints[END_REF] minimize the makespan of a flow shop under peak power consumption constraints. GHG emission management is another major challenge at the operational level. [START_REF] Foumani | The impact of various carbon reduction policies on green flowshop scheduling[END_REF] provide a comprehensive study on the different carbon taxes policies applied to a flow shop. [START_REF] Zhang | Energy-conscious flow shop scheduling under time-of-use electricity tariffs[END_REF] develop a time-indexed formulation for a flow shop under TOU tariffs and carbon emission periods with a trade-off between the low-carbon emission period and the TOU onpeak hours. In their work, peak demands are handled by natural gas and base energy is provided by coal-based sources, which emit more GHG.

Few studies have been done for the OAS under TOU tariffs, and even fewer have been focused on GHG emissions. The work of [START_REF] Chen | Order acceptance and scheduling problem with carbon emission reduction and electricity[END_REF] is the first of its kind, while proposing a benchmark and an exact solving approach with a disjunctive formulation. On this basis and upon their study, Section 3 provides two time-indexed MILP for the OAS problem with release dates under energy aspects. In Section 4 and 5, a comparative analysis based on the number of feasible and optimal solutions found and the time spent to solve instances is carried out. In addition, the characteristics of the models in terms of number of constraints and variables are discussed.

PROBLEM FORMULATION AND SOLUTION APPROACH

We investigate an OAS problem on a single nonpreemptive machine with release dates. Each order 𝑗 is characterized by its processing time 𝑝 𝑗 , release date 𝑟 𝑗 , due date 𝑑 𝑗 , deadline 𝑑 ̅ 𝑗 , revenue 𝑒 𝑗 , tardiness penalty 𝑤 𝑗 and power consumption Ω 𝑗 .

Figure 1: Profit calculation of an order 𝑗 taken from [START_REF] Chen | Order acceptance and scheduling problem with carbon emission reduction and electricity[END_REF] Figure 1 presents the profit calculation of an order 𝑗 according to its deadline, due date and tardiness penalties. An order 𝑗 is accepted only if it can be finished before its deadline 𝑑 ̅ 𝑗 ; besides , it earns the totality of its revenue 𝑒 𝑗 until its due date 𝑑 𝑗 , since the tardiness penalties 𝑤 𝑗 are applied after 𝑑 𝑗 .

In the initial work of [START_REF] Chen | Order acceptance and scheduling problem with carbon emission reduction and electricity[END_REF], the horizon is divided into 𝑚 and ℎ intervals of electricity TOU and carbon emission respectively. Each electricity TOU interval 𝑘 = 1, … , 𝑚 is characterized by an electricity cost 𝐸𝐶 𝑘 and a starting time 𝑏 𝑘 and each carbon emission interval 𝑙 = 1, … , ℎ is defined by its starting time 𝑔 𝑙 and an amount 𝑞 𝑙 of emitted carbon per kg, applying a 𝑇𝑎𝑥 for each emitted kg of CO2. For the sake of simplicity, in our formulations the time horizon is split into 𝑇 periods determined in equation ( 1), where each period is characterized by its electricity TOU cost and the amount of emitted CO2.

𝑇 = max 𝑗=1,…,𝑛 𝑑 ̅ 𝑗 + 1 (1)
In this problem, idle times are considered but their energy consumption is neglected. The objective is to maximize the sum of the profit of the orders minus the environmental costs during processing time.

Two time-indexed MILP models are developed for this profit and time-driven problem. These formulations rely on the discretization of time, i.e. time is divided into unitary slots 𝑡 = 0, … , 𝑇.

The energy cost 𝑐 𝑗𝑡 is precomputed for each order 𝑗 = 1, … , 𝑛 and at each time processing time 𝑡 = 𝑟 𝑗 , … , 𝑑 ̅ 𝑗 :

𝑐 𝑗𝑡 = Ω 𝑗 60 (∑ 𝐸𝐶 𝑘 𝕝𝑡≥𝑏 𝑘-1 𝑡<𝑏 𝑘 𝑚 𝑘=1 + 𝑇𝑎𝑥 ∑ 𝑞 𝑙 𝕝𝑡≥𝑔 𝑙-1 𝑡<𝑔 𝑙 ℎ 𝑙=1 ) (2)
In equation ( 2

Pulse formulation

The first time-indexed model is known as the pulse formulation, where the binary decision variables 𝑥 𝑗𝑡 = 1 indicates that order 𝑗 = 1, … , 𝑛 starts at time 𝑡 = 1, . . , 𝑇 or not (𝑥 𝑗𝑡 = 0). Note that the possible periods for the starting time of an order 𝑗 are 𝑡 = 𝑟 𝑗 , … , 𝑑 ̅ 𝑗 -𝑝 𝑗 + 1.

The profit 𝑓 𝑗𝑡 of an order 𝑗 = 1, … , 𝑛 at time 𝑡 = 𝑟 𝑗 + 𝑝 𝑗 , … , 𝑑 ̅ 𝑗 is precomputed. Equation (3) corresponds to this calculation, which is the revenue minus the possible tardiness penalties when 𝑡 > 𝑑 𝑗 .

𝑓 𝑗𝑡 = 𝑒 𝑗 -𝑤 𝑗 max{0; 𝑡 -𝑑 𝑗 }

(3)

The corresponding MILP is written as follows.

maximize ∑ ∑ 𝑥 𝑗𝑡 (𝑓 𝑗𝑡 -(∑ 𝑐 𝑗𝑡 ′ 𝑡+𝑝 𝑗 -1 𝑡 ′ =𝑡
))

𝑑 ̅ 𝑗 -𝑝 𝑗 +1 𝑡=𝑟 𝑗 𝑛 𝑗=1

(4)

∑ 𝑥 𝑗𝑡 𝑛 𝑗=1 ≤ 1, ∀𝑡 = 0, … , 𝑇 (5)
∑ 𝑥 𝑗𝑡 𝑑 ̅ 𝑗 -𝑝 𝑗 +1 𝑡=𝑟 𝑗 ≤ 1, ∀𝑗 = 1, … , 𝑛 (6)
∑ 𝑥 𝑗𝑡 𝑟 𝑗 -1 𝑡=0 = 0, ∀𝑗 = 1, … , 𝑛 (7)
∑ 𝑥 𝑗𝑡 𝑇 𝑡=(𝑑 ̅ 𝑗 -𝑝 𝑗 +1)+1 = 0, ∀𝑗 = 1, … , 𝑛 (8)
𝑥 𝑖𝑡 + ∑ 𝑥 𝑗𝑡 ′ 𝑡+𝑝 𝑖 -1 𝑡 ′ =𝑡 𝑟 𝑗 +1≤𝑡 ′ ≤𝑑 ̅ 𝑗 ≤ 1, ∀𝑖, 𝑗 = 1, … , 𝑛; 𝑖 ≠ 𝑗, ∀𝑡 = 𝑟 𝑖 , … , 𝑑 ̅ 𝑖 -𝑝 𝑖 + 1 (9)
The objective ( 4) is the maximization of the sum of the total profit of each order, that is the profit 𝑓 𝑗𝑡 including the tardiness penalties and environmental cost during the processing time, given by the sum of 𝑐 𝑗𝑡 ′ from the starting time 𝑡 until completion 𝑡 + 𝑝 𝑗 -1. Constraints (5) specify that at each time 𝑡, the machine can start only one job. Constraints (6) restrict the starting time of each order to the interval defined from its release date to its deadline.

In the same manner, constraints ( 7) and ( 8) prevent each order to be processed before its release date and after its deadline. Constraints (9) prevent any order 𝑗 to overlap in the interval [𝑡, 𝑡 + 𝑝 𝑖 -1] when an order 𝑖 starts at time 𝑡 = 𝑟 𝑖 , … , 𝑑 ̅ 𝑖 -𝑝 𝑖 + 1.

Table 1 presents the optimal solution of an example with 𝑛 = 4 orders with their processing times 𝑡 𝑗 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 4 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 Table 1: Optimal solution represented by the values of the decision variables 𝑥, sequence is 4-2-3 and order 1 is rejected.

For this example, there are 180 constraints and 80 variables. The optimal solution has been found in 0.01 s.

Figure 2: Gantt chart representation for the example.

On-off formulation

The second model presented is the on-off formulation.

Each binary decision variable 𝑥 𝑗𝑡 = 1 indicates whether the order 𝑗 is processed at time 𝑡 = 𝑟 𝑗 , … , 𝑑 ̅ 𝑗 , or not 𝑥 𝑗𝑡 = 0. In addition, the binary decision variable 𝑎 𝑗 = 1 represents whether the order 𝑗 is accepted.

maximize ∑ (𝑎 𝑗 𝑒 𝑗 -𝑤 𝑗 max 𝑡=𝑑 𝑗 ,…𝑑 ̅ 𝑗 {𝑡 -𝑑 𝑗 }𝑥 𝑗𝑡 - 𝑛 𝑗=1 ∑ 𝑐 𝑗𝑡 𝑥 𝑗𝑡 𝑑 ̅ 𝑗 𝑡=𝑟 𝑗 ) (10) ∑ 𝑥 𝑗𝑡 𝑛 𝑗=1 ≤ 1, ∀𝑡 = 0, … , 𝑇 (11)
∑ 𝑥 𝑗𝑡 𝑟 𝑗 -1 𝑡=0 = 0, ∀𝑗 = 1, … , 𝑛 (12)
∑ 𝑥 𝑗𝑡 𝑇 𝑡= 𝑑 ̅ 𝑗 +1 = 0, ∀𝑗 = 1, … , 𝑛 (13)
𝑝 𝑗 𝑎 𝑗 = ∑ 𝑥 𝑗𝑡 𝑑 ̅ 𝑗 𝑡=𝑟 𝑗 , ∀𝑗 = 1, … , 𝑛 (14)
∑ 𝑥 𝑗𝑡 ′ 𝑡-𝑝 𝑗 𝑡 ′ =𝑟 𝑗 + ∑ 𝑥 𝑗𝑡 ′ 𝑑 ̅ 𝑗 𝑡 ′ =𝑡+𝑝 𝑗 ≤ (1 -𝑥 𝑗𝑡 )𝑝 𝑗 , ∀𝑗 = 1, … , 𝑛, ∀𝑡 = 𝑟 𝑗 , … , 𝑑 ̅ 𝑗 ( 15
)
The objective function ( 10) is the maximization of the total profit of the accepted orders, i.e. their revenues 𝑎 𝑗 𝑒 𝑗 minus their possible tardiness penalties 𝑤 𝑗 by retrieving the instant 𝑡 = 𝑑 𝑗 , … , 𝑑 ̅ 𝑗 when 𝑥 𝑗𝑡 = 1 (completion time) and the environmental costs during processing time. Constraints (11) state that at each time the machine is either doing nothing or processing an order. Constraints (12) and ( 13) ensure that each order cannot be processed before its release date and after its deadline. Constraints (14) impose that each accepted order must be processed during the totality of its processing time. Constraints (15) guarantee non-preemption by forcing the continuity of the decision variables.

MOSIM '20 -November 12-14, 2020 -Agadir -Morocco Table 2 presents the same solution as in 3.1 with the perspective of the on-off model.

𝑡 𝑗

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 4 0 0 1 1 1 1 0 0 0 0 0 0 0 0 Table 2: Optimal solution represented by the values of the decision variables 𝑥, sequence is 4-2-3.

For this solution and this formulation, there are 97 constraints and 111 variables. The optimal solution has been found in 0.03s.

RESOLUTION METHOD

The considered benchmark derives from the work of [START_REF] Chen | Order acceptance and scheduling problem with carbon emission reduction and electricity[END_REF]. Similar TOU tariffs and carbon emission periods are used. Instances with various number of orders 𝑛 = 5, 10, 15, 20, 50 and with different characteristics are tested in order to have a diverse set of instances. Two coefficients are used to generate the parameters of an instance: the tardiness factor 𝜏 = 0.1, 0.3, 0.5 and the due date ranges 𝑅 = 0.1, 0.5, 0.9. This first comparative study involves 45 instances.

Processing time, release date and revenue are generated with an uniform distribution: 𝑝 𝑗 , 𝑒 𝑗 ~𝒰(1,20) and 𝑟 𝑗 ~𝒰(0, 𝜏𝑝 𝑇 ) with 𝑝 𝑇 = ∑ 𝑝 𝑗 𝑛 𝑗=1

. Due dates, deadlines, tardiness penalties and power consumption are computed from the values of these generated parameters:

𝑑 𝑗 = 𝑟 𝑗 + max{𝑠𝑙𝑎𝑐𝑘, 𝑝 𝑗 } (16) 𝑠𝑙𝑎𝑐𝑘 ~ 𝒰 (𝑝 𝑇 (1 -𝜏 - 𝑅 2 ) , 𝑝 𝑇 (1 -𝜏 + 𝑅 2 )) (17)
𝑑 ̅ 𝑗 = 𝑑 𝑗 + 𝑅𝑝 𝑗 (18)

𝑤 𝑗 = 𝑒 𝑗 𝑑 𝑗 -𝑑 𝑗 (19) Ω 𝑗 ~𝒰(1, 𝑒 𝑗 ) × 1 2 (20)
This ensures coherent values for these parameters. For instance, the deadline must be greater or equal to the due date.

The tested models are the pulse formulation, the on/off formulation and the disjunctive model of [START_REF] Chen | Order acceptance and scheduling problem with carbon emission reduction and electricity[END_REF]. For the sake of comparability, each model has been implemented and solved in a commercial solver (IBM CPLEX Optimization Studio v12.9) on a desktop computer with processor Intel i5 2GHz CPU with 4GB RAM. Solving time is limited to 3600 seconds.

COMPUTATIONAL RESULTS

Table 3 resumes the benchmark results for each model and their average performances for 9 instances of same size and different values for 𝜏 and 𝑅. Average solving time (𝑐𝑝𝑢), average gap (𝑔𝑎𝑝) and the number of feasible (#𝑓𝑒𝑎) and optimal (#𝑜𝑝𝑡) solutions found are reported.

In our tests, the gap is retrieved from CPLEX relative MILP gap, which represent the gap between the best bound and the best integral solution found by the solver. A summary of the performances of all formulation on the 45 instances is also provided.

The results clearly point out that time-indexed formulations outperform a standard disjunctive model on average. The pulse formulation find 38 optimal solutions among 44 in 711 seconds on average, which represent a success rate of 86%. In contrast, the disjunctive formulation takes on average more than twice the time to find half of the optimal solutions. The on/off formulation provides similar results to the pulse model, achieving a rate of 84% of optimal solutions found in 721 seconds on average. This model seems to be more performant than the pulse formulation for small to medium instances (𝑛 = 5, 10, 15, 20), finding all the optimal solutions in less time. However, its average gap is less tight than the pulse formulation.

The size of time-indexed formulations (constraints and variables) is their main weakness. Indeed, both timeindexed formulations cannot find all the feasible solutions of the benchmark due to either out-of-memory issues or low quality gap. At worst, the on-off formulation has 𝒪(𝑛𝑇) + 𝒪(𝑛) variables and 𝒪(𝑛𝑇) constraints, whereas the pulse formulation has 𝒪(𝑛𝑇) variables and 𝒪(𝑛 2 𝑇) constraints. The disjunctive formulation has 𝒪(𝑛 2 ) + 𝒪(𝑛𝑚) + 𝒪(𝑛ℎ) variables and constraints.

As seen in Section 3, the pulse formulation binary variables contain information for both acceptance and the instant where the order starts, whereas the on/off formulation differs in semantic. A binary variable 𝑥 𝑗𝑡 in the pulse formulation corresponds to 𝑝 𝑗 + 1 variables in the on/off formulation. The pulse formulation effectiveness on average compared to the on/off model may reside in the use of less binary variables and thus imply a lower number of branching. Nevertheless, the number of constraints in the on/off model provides an advantage for small instances.

CONCLUSION AND PERSPECTIVES

In this paper, we presented two time-indexed formulations for the OAS with release dates under energy aspects.

Our proposed formulations are more performant on average than the disjunctive formulation described in [START_REF] Chen | Order acceptance and scheduling problem with carbon emission reduction and electricity[END_REF]. This can be explained by their LP-relaxation which provide good bound for medium instances, according to [START_REF] Van Den Akker | Time-Indexed Formulations for Machine Scheduling Problems: Column Generation[END_REF]. Moreover, time-indexed models seem to be the most efficient formulations '20 -November 12-14, 2020 -Agadir -Morocco for this problem since the objective is the maximization of a time-driven profit comprising time-related penalties and time-varying environmental costs. However, these formulations are limited by their spatial complexity. Our future work is focused on the development of dedicated cuts or other exact approaches that have the potential to improve the performances for large instances. Finally, time-indexed formulations for an extension of the proposed problem with setup-dependent sequence are under development.

MOSIM

summary

  ), 𝕝 𝑥 denotes the indicator function, which takes value 1 if condition 𝑥 holds and 0 otherwise. For each order 𝑗 at each time 𝑡, the power consumptionexpressed into minutesis multiplied by the cost of the respective TOU interval 𝑘 = 1, … , 𝑚 of period 𝑡 ∈ [𝑏 𝑘-1 , 𝑏 𝑘 [ and the taxed CO2 emission interval 𝑙 = 1, … , ℎ of period 𝑡 ∈ [𝑔 𝑙-1 , 𝑔 𝑙[ . MOSIM'20 -November 12-14, 2020 -Agadir -Morocco

  𝑝 = [5,3,2,4], release dates 𝑟 = [1,2,1,1], due dates 𝑑 =[6,5,12,7], deadlines 𝑑 ̅ =[9,10,14,12], revenues 𝑒 = [10,10,6,10], power consumption Ω = [1,2,1,1] and weight penalties 𝑤 = [2,1,3,2]. Finally, the starting times of TOU and carbon emission intervals 𝑏 = 𝑔 = [0,5,8], the electricity price 𝐸𝐶 = [2,10,2] and the amount of CO2 emitted 𝑞 =[START_REF] Liao | Multi-objective optimization of single machine scheduling with energy consumption constraints[END_REF]1,[START_REF] Liao | Multi-objective optimization of single machine scheduling with energy consumption constraints[END_REF] are defined.

Table 3 :

 3 Models performances in terms of solving time, gap to the best-bound, feasible and optimal solutions (*out-of-memory status or found infinite gap).

	45	24	1693 2.5	44	38	711 0.05 44	37	721	1.2

ACKNOWLEDGMENTS

This research was funded by the Grand-Est region and the Aube department in France.


Xuan Wu
email: xuan.wu@ensam.eu

Virginie Goepp
email: virginie.goepp@insa-strasbourg.fr

Ali Siadat
email: ali.siadat@ensam.eu

FUNCTION ANALYSIS OF MANUFACTURING EXECUTION SYSTEMS (MES) IN CYBER PHYSICAL PRODUCTION SYSTEMS (CPPS)

Keywords: Cyber physical systems, Cyber physical production systems, Manufacturing execution systems, Functional modelling (IDEF0), 5C architecture, ANSI/ISA-95

published or not. The documents may come

INTRODUCTION

In the era of Industry 4.0, Cyber Physical Systems (CPSs) play an essential role because they make the fusion of the physical and the virtual worlds [START_REF] Kagermann | Securing the future of German manufacturing industry: Recommendations for implementing the strategic initiative industrie 4.0[END_REF]. The application of CPSs in production environments leads to the development of Cyber Physical Production Systems (CPPSs). The first detailed description of a CPPS was given by [START_REF] Monostori | Cyber-physical systems in manufacturing[END_REF], which has been widely cited and broadly accepted in recent years. As a complement to Monostori's description, Cardin (2019) added several missing points, including knowledge management, decision making and adaptability, and adapted it as following (the complementary part of Cardin's description has been highlighted in bold): "CPPSs are systems of systems of autonomous and cooperative elements connecting with each other in situation dependent ways, on and across all levels of production, from processes through machines up to production and logistics networks, enhancing decisionmaking processes in real-time, response to unforeseen conditions and evolution along time". [START_REF] Monostori | Cyber-physical systems in manufacturing[END_REF] claimed that CPPSs partly break with the traditional automation pyramid. Traditionally, the Manufacturing Execution Systems (MES) was positioned as a middle layer to bridge the shop-floor with the Enterprise Resource Planning (ERP) system (left side of Figure 1). Today, the automation pyramid will be transformed into a more decentralized structure (right side of Figure 1). However, this transformation leads to vague on the role of MES in CPPSs: Does this decentralized structure of CPPSs make MES superfluous? If not, what is the role of MES in CPPSs? And how does MES need to change to fulfill this role? Some scholars attempted to answer these questions. For example, [START_REF] Rossit | Industry 4.0: Smart Scheduling[END_REF] claimed that current MES, which take care of scheduling and dispatching work orders, will be absorbed by CPPSs. This minimizes response times, increase the flexibility of the entire system, which is useful for decision making. However, [START_REF] Arica | Status and future of manufacturing execution systems[END_REF] and [START_REF] Mantravadi | An Overview of Next-generation Manufacturing Execution Systems: How important is MES for Industry 4.0?[END_REF] outlined that with data collection, analysis, and communication functions across the value chains, MES should play a central role in Industry 4.0 manufacturing systems. ZVEI, (2017) working group also indicated that the significance of MES in Industry 4.0 will increase and not decrease. But the focus will shift away from simple execution management towards comprehensive coverage of all Manufacturing Operation Management (MOM) activities and away from the control of production towards optimization including the incorporation of events. Almada-Lobo, (2016) pointed out that MES is required to change to cope with the challenges created by Industry 4.0, especially considering the following four aspects: decentralization, vertical integration, connectivity and cloud computing and advanced analysis. We find that, in current academic research, there still remains a dispute on the role of MES in CPPSs, ranging from "no role because the hierarchical structure will be dissolved" and "fewer and fewer role because MES functions will be cannibalized by enhanced CPPSs and ERP functions" right up to "central role because MES will be the decision-making MOSIM '20 -November 12-14, 2020 -Agadir -Morocco center of an organization and create the optimal value chain". The existing studies are either based on literature reviews or empirical methods, but no one builds functional models to elaborate and compare the functions between CPPSs and MES. Therefore, this paper aims to shed light on the role of MES in CPPSs, by proposing functional models of CPPSs and MES in the form of IDEF0 diagrams and thus comparing their functions. The rest of the paper is organized as follows. In Section 2, functional models of CPPSs and MES are proposed using IDEF0 diagrams. In Section 3, the data flows between the functional models of CPPSs and MES are described and then the functions of the two are compared. Section 4 concludes this paper and indicates future work.

FUNCTIONAL MODELS OF CPPS AND MES

IDEF0 (Integration DEFinition of function modelling) offers a functional modelling methodology for activity and information-flow analysis, focusing on the hierarchical decomposition of activities and the interaction flows between activities (IDEF0, 1993). In this section, using IDEF0 diagrams, the functional model of a CPPS based on the 5C architecture proposed in [START_REF] Lee | A cyberphysical systems architecture for Industry 4.0-based manufacturing systems[END_REF] and the functional model of an MES based on the ISA-95 standard (ISA-95 2005) are worked out. Level IV (cognition level) generates a thorough insight of the whole system, which supports expert users to make decisions. At Level V (configuration level), CPPSs can self-configure and self-adapt to return to normal behaviors. For each level, we can extract a major activity, which is "acquire raw data from the physical space", "convert raw data to meaningful data", "build the cyber space", "make decisions", and "self-configure and self-adapt", respectively.

Functional models of CPPSs

5C architecture

Functional model of CPPSs

According to the description of the 5C architecture, each level has a major activity. Therefore, the context diagram A0 "Implement a CPPS" can be decomposed into these five sub-activities, A1-A5, as shown in Figure 3, which shows the activities involved in the manufacturing process along with the information flows required to support these activities. A detailed description of the function modules and their relationship is given as follows.

A1-Acquire raw data from the physical space

This box describes the main activity in the Level I of the 5C architecture. The input is "operational commands", which is transformed into "raw data". The raw data refers to measurements data from plant instrumentation. The activity uses the mechanisms "data sources (such as sensors and controllers)", "communication protocols (such as MTConnect)" and "data acquisition methods", with the controls "production rules", "production capability" and "production routing".

A2-Convert raw data to meaningful data

This box describes the main activity in the Level II of the 5C architecture. The input that is linked to the prior output "raw data" of the activity A1 is transformed into "meaningful data". The meaningful data refers to the data at a higher semantic level, such as movement events of resources and statuses data of resources. Data analysis technologies as the mechanism are needed for getting meaningful data. In addition to having the same controls as the activity A1, another control is the data format.

MOSIM'20 -November 12-14, 2020 -Agadir -Morocco A3-Build the cyber space This box describes the main activity in the Level III of the 5C architecture. It has one input which is linked to the prior output "meaningful data" of the activity A2, one output "simulation data of status and behaviors of CPPSs", and the same controls as the activity A1. Simulation software and information systems are needed to support the activity and its function.

A4-Make decisions

This box describes the main activity in the Level IV of the 5C architecture. It has one input which is linked to the prior output "simulation data" of the activity A3, and one output "corrective and preventive decisions". Simulation software, information systems, proper presentation tools and expert users are needed to support the activity and this function. In addition to having the same controls as the activity A1, another control, the risk criterion, is needed to support the correct decisions to be taken.

A5-Self-configure and self-adapt

This box describes the main activity in the Level V of the 5C architecture. The input that is linked to the prior output "decisions" of the activity A5 is transformed into "adjusted production activities". It has the same controls as the activity A1. Physical entities in shop floor are needed to support the activity and its function.

Functional models of MES

ANSI/ISA-95 standard

The part 3 of the standard ISA-95 (2005) gives a detailed overview of all the activities and data flows of MOM in the area of production, inventory, maintenance and quality operational management. MOM could be regarded as the upper limited range of functions involved in MES applications. Generally speaking, the main focus of current MES applications is only on the production operations management. However, future MES functions will have to be developed towards MOM. Therefore, in this paper, MES refers to MOM, which covers all areas of operations management. One of the main advantages of ISA95 is the definition of the main functions and data flows within MOM. Therefore, for this study, MES function's classification is based on the ISA-95.

Functional model of MES

The ISA-95 described the generic activities in MOM, as shown in Figure 4. Each of the four areas of operations management (production, inventory, maintenance and quality) is composed of these eight activities. According to this, the context diagram B0 "Implement an MES" was decomposed into eight sub-activities, B1-B8, as depicted in Figure 5. Activity B1 (Manage definitions of work) manages the work information and develop work rules. Activity B2 (Manage resource) manages the information about the resources including machines, tools, labor, materials and MOSIM'20 -November 12-14, 2020 -Agadir -Morocco energy. Activity B3 (Detailed schedule) makes the schedule and determines the optimal use of resources to meet the schedule requirements. Activity B4 (Dispatch work) dispatches work lists to equipment and personnel. Activity B5 (Manage executions of work) directs the performance of work, as specified by the dispatch list. Activity B6 (Track work) summarizes and reports information about the personnel and equipment used, the material consumed, material produced and other relevant data such as costs and performance analysis results. Activity B7 (Collect data) gathers, compiles and manages data for specific work processes. Activity B8 (Analyze performance) analyzes and reports performance information of production systems.

The data flows between these activities are extracted from part 3 of the ISA-95. For example, one of the tasks in product definition management, a subset of "B1manage definitions of work", is extracted from the ISA-95: providing product production rules to personnel or other activities. From this task, we can conclude that the production rule is the output of B1. Like this, the data flows between these activities can be found in the ISA-95, so we will not discuss them here.

FUNCTION ANALYSIS OF MES IN CPPS

Based on the established functional model of CPPSs and MES, this section describes the data flows between CPPSs and MES, which represents the functional interfaces between them. Then based on the data flow, this section analyses which functions of MES are redundant in CPPSs and which functions are central.

There are two-way data flows in MES. On the one hand, MES deals with the top-down data flow: the requirements provided by the organizational level must be transformed into detailed operational commands to the production areas. On the other hand, it manages the bottom-up data flow. MES collects shop floor data, analyzes it, and extracts useful information to provide key performance indicators (KPIs) to the company for making commercial decisions and improving production performance. Therefore, the data flows between CPPSs and MES are described as follows.

 Data flow 1: "equipment and process rules" from the output of the Activity B1  mechanisms of the Activity A1 to A5.

Equipment and process rules define the specific production/inventory/maintenance/quality test instructions that are executed in the shop floor based on the specific assigned tasks, such as programs for CNC machines.  Data flow 2: "operational commands" from the output of the Activity B5  the input of the Activity A1.

Operational commands define the commands to start or complete elements of a work order in the shop floor, such as procedures for setting up machines or maintenance of machines.  Data flow 3: "operational responses" from the output of Activity A3  the input of the Activity B5.

Operational responses define information received from the shop floor in response to commands, such as the completion or status of elements of work orders.  Data flow 4: "equipment and process specific data, resource data, maintenance history, quality test data, inventory data " from the output of Activity A3  the input of the Activity B7.

Equipment and process specific data defines the data about the process being performed and the resources involved.  Data flow 5: "Production, maintenance, quality and inventory related KPIs" from the output of Activity B8

 the input of Activity A4

The manufacturing operations management areas have different sets of performance indicators, which are used together to make decisions and to monitor the realization of enterprise business objectives.  Data flow 6: "current resource availability" from the output of Activity A3  the input of Activity B2.

Resources include machines, tools, labor (with specific skill sets), materials, and energy.

 Data flow 7: "resource capability" from the output of Activity B2  the input of A3

Resource capability defines the committed available, or unattainable of resources. The information is based on the current statuses and future needs (as identified in Activity B3).  Data flow 8: "production and movement events" from the output of Activity A3  the input of the Activity B6.

Production and movement events include the movements of the path of materials and updates to lot locations.

For the data flow 1, 3, 4, 6, 7 and 8, MES mainly serves as an information platform that contains all information in the production, inventory, maintenance and quality test operations management and we find that the functions of MES are very similar to those of digital twin in terms of work definitions management, resource management, work track and data collection. For example, in terms of the MES's function "data collection", digital twins store all data from the field and other information systems in real-time, as in the MES. Digital twin refers to the virtualization of the physical resources in the cyber world and its main purpose is to facilitate the decision-making process through real-time simulation [START_REF] Tao | Digital twin-driven product design, manufacturing and service with big data[END_REF]. Therefore, future MES can be used as a digital twin due to its ability to provide digital images and monitor manufacturing processes in real-time.

For the data flow 5, MES analyzes and reports performance information to expert users for decision-making.

In the Industry 4.0 context, the function of performance analysis in MES becomes more complex and important. With the development of information and communication technology (ICT), the amount of data collected from the shop floor increases tremendously. MES should always be required to present the most relevant data to the humans for decision making.

For the data flow 2, MES sends the operational commands to the physical space of CPPSs. The function of executions management is equally important because it is not only necessary to have data coming from the shop floor to MES, but also to be able to act from the MES to the shop floor, by offering operational commands on the production control system. In the CPPSs context, it should be possible in the future to decentralize operational decisions and delegate a part of decisions made previously in MES to lower levels of production systems, such as smart machines and smart products. This reduces the number of tasks of this function in terms of the operational level decisions. Since MES covers all the areas of MOM, it can collect all data related to production, quality, inventory and maintenance. This endows MES with the holistic insights for the manufacturing processes and, which can make commands before any disruptions.

CONCLUSION AND FUTURE WORK

The paper studies the role of MES in CPPSs by analyzing the functional models of CPPSs and MES. MES's functions are still required for CPPSs, but they need to MOSIM'20 -November 12-14, 2020 -Agadir -Morocco be presented in new and different forms, with some being incorporated into the CPPSs' functions, some expanding their horizons and some decentralizing to the local intelligent entities. For example, the MES's function "detailed schedule" will only have high level tasks, making rough planning of manufacturing orders, while flexible planning of partial orders will be made by local entities of CPPSs through taking into account local resource availabilities in real-time. Correspondingly, MES software will move from monolithic software systems towards functional modules that will be implemented in a highly distributed way. Therefore, the next-generation MES requires greater flexibility, dynamism and improved functionality than today's MES.

In the future work, validation in a real platform will be worked out to evaluate the quality of the proposed model. Besides, we will implement a "dynamic link" between CPPSs and MES from a functional perspective. In our view, the functional level is the best suited for this implementation as it enables to free of the underlying technical solutions. Indeed, the MES and CPPS interrelated functions cannot be engineered independently and then simply combined together. The idea is rather to coordinate the design of both systems to explore hidden functions and ensure that their interactions are mutually supportive. This "dynamic link" is twofold: one the one hand, if the configuration of CPPSs (e.g. the layout of CPPSs) changes, the organization and business processes will be changed automatically. One the other hand, if new functionalities of EISs are needed, CPPSs can configure themselves automatically to meet these new requirements. The digital twin could be a solution to implement such dynamic link as the change of state of the physical object will have impacts on the state of the cyber object automatically and vice versa in digital twins. [START_REF] Cimino | Review of digital twin applications in manufacturing[END_REF] pointed out the future possibility of integrating the digital twin of the field devices with the MES in order that decision making could be made in a fully reactive way through the digital twin to the MES.

Figure 1 :

 1 Figure 1: Decomposition of the automation hierarchy with distributed services (VDI/ VDE, 2013)

Figure 2 :

 2 Figure 2: 5C architecture for the implementation of CPSs (Lee et al., 2015)

Figure 3 :

 3 Figure 3: Activity A0 -Implement a CPPS.

Figure 4 :

 4 Figure 4: Generic activity model of manufacturing operations management (ISA-95, 2005)

ACKNOWLEDGMENTS

The paper of the first author is financially supported by China Scholarship Council under Grant 201706020154.