Shiling Feng

Xuanbing Qiu

Guqing Guo

Enhua Zhang

Qiusheng He

Xiaohu He

Weiguang Ma

Christa Fittschen

Chuanliang Li
email: li_chuanliang@126.com

A palm-sized laser spectrometer with high robustness and sensitivity for trace gas detection using a novel double-layer toroidal cell

Keywords: palm-sized laser spectrometer, double-layer toroidal cell, near infrared spectroscopy, calibration-free scanned WMS, vibration and temperature environment tests

A palm-sized laser spectrometer has been developed for detecting trace gases based on tunable diode laser absorption spectroscopy (TDLAS) in combination with a novel double-layer toroidal cell. With the benefit of a home-made electronic system and compact optical design, the physical dimension of the sensor is minimized into 24×15×16 cm 3 . A toroidal absorption cell, with 84 reflections in 2 layers for an effective optical path length of 8.35 m, was used to enhance the absorption signals of gaseous species. A home-made electronic system was designed for implementing distributed feedback (DFB) diode laser controller, analog lock-in amplifier, data acquisition and communication. Calibration-free scanned wavelength modulation spectroscopy was employed to determine the concentration of the gas and reduce the random fluctuations from electronical noise and mechanical vibration. The measurement of CH4 in ambient air was demonstrated using a DFB laser at 1.653 μm. The rise time and fall time for renewing the gas mixture are approximately 16 s and 14 s, respectively. Vibration and temperature tests have been carried out for verifying the performance of the spectrometer, and standard deviations of 0.38 ppm and 0.11 ppm for 20 ppm CH4 at different vibration frequencies and temperatures, respectively, have been determined. According to the Allan deviation analysis, the minimum detection limit for CH4 can reach 22 ppb at an integration time of 57.8 s. Continuous measurement of atmospheric CH4 for two days validated the feasibility and robustness of our laser spectrometer, providing a promising laser spectral sensor for deploying in unmanned aerial vehicles or mobile robots.

Introduction

Trace gas detection plays an important role in multiple application fields, such as analyzing human breath [START_REF] Heinrich | Infrared laser-spectroscopic analysis of 14 NO and 15 NO in human breath[END_REF] , monitoring pollutants 2 or greenhouse gases [START_REF] Wang | Applications of a tunable diode laser absorption spectrometer in monitoring greenhouse gases[END_REF] and measuring gases of interest in industrial process [START_REF] Hirschmann | Trace gas detection of benzene, toluene, p-, m-and o-xylene with a compact measurement system using cantilever enhanced photoacoustic spectroscopy and optical parametric oscillator[END_REF] . A variety of conventional techniques are successfully used for trace gas detection, such as gas chromatography [START_REF] Loftfield | Automated Gas Chromatographic System for Rapid Analysis of the Atmospheric Trace Gases Methane, Carbon Dioxide, and Nitrous Oxide[END_REF] , mass spectroscopy [START_REF] Maity | Formation of Ketene (H2CCO) in Interstellar Analogous Methane (CH4)-Carbon Monoxide (CO) ICES: A Combined Ftir and Reflectron Time-of-Flight Mass Spectroscopic Study[END_REF] and Fourier transform spectroscopy [START_REF] Tanaka | Aircraft measurements of carbon dioxide and methane for the calibration of ground-based high-resolution Fourier Transform Spectrometers and a comparison to GOSAT data measured over Tsukuba and Moshiri[END_REF] . Although these methods have high sensitivity and selectivity, they are costly, time-consuming and bulky. Therefore, these approaches are not suitable for a real-time, in-situ and mobile measurement of trace gases. Tunable diode laser absorption spectroscopy (TDLAS) has been widely employed in detecting trace gases due to its high selectivity, sensitivity and rapid response [START_REF] Pogany | Development and metrological characterization of a tunable diode laser absorption spectroscopy (TDLAS) spectrometer for simultaneous absolute measurement of carbon dioxide and water vapor[END_REF][START_REF] Zellweger | Evaluation of new laser spectrometer techniques for in-situ carbon monoxide measurements[END_REF][START_REF] Dong | Compact TDLAS based sensor design using interband cascade lasers for mid-IR trace gas sensing[END_REF][START_REF] Wienhold | TRISTAR-atracer in situ TDLAS for atmospheric research[END_REF] . In recent years, the trend of compact laser spectrometers has been the focus of TDLAS for fulfilling the embedding requirements for unmanned aerial vehicles [START_REF] Yang | Natural Gas Fugitive Leak Detection Using an Unmanned Aerial Vehicle: Measurement System Description and Mass Balance Approach[END_REF] or mobile robots [START_REF] Bennetts | Towards Real-World Gas Distribution Mapping and Leak Localization Using a Mobile Robot with 3D and Remote Gas Sensing Capabilities[END_REF] . Due to restriction of physical dimension, the effective optical path length is commonly shortened, which decreases the detection sensitivity of the sensor. Therefore, a multipass cell (MPC) is a better choice for compact laser spectrometer. White cells [START_REF] White | Long Optical Paths of Large Aperture[END_REF][START_REF] Pilston | A Long Path Gas Absorption Cell[END_REF] , Herriott cells [16][START_REF] Guo | A portable sensor for in-situ measurement of ammonia based on near-infrared laser absorption spectroscopy[END_REF][START_REF] Li | High-speed multi-pass tunable diode laser absorption spectrometer based on frequency-modulation spectroscopy[END_REF] and Chernin absorption cells [START_REF] Chernin | Optical multipass matrix systems[END_REF] are three typical MPC used in TDLAS. These absorption cells consist of two or more individual mirrors. However, they are susceptible to mechanical vibration and temperature variation. Thus, a compact, stable and welladapted absorption cell is desirable for a TDLAS senor to be employed in artificial intelligence inspection vehicles. Recently, a compact toroidal MPC with low-volume monolithic ring [START_REF] Tuzson | Compact multipass optical cell for laser spectroscopy[END_REF][START_REF] Jouy | Mid-infrared spectroscopy for gases and liquids based on quantum cascade technologies[END_REF] was presented and we also proposed a toroidal cell with multi-layer patterns [START_REF] Chang | Implementation of the toroidal absorption cell with multi-layer patterns by a single ring surface[END_REF] .

Wavelength modulation spectroscopy (WMS) [START_REF] Guo | A portable sensor for in-situ measurement of ammonia based on near-infrared laser absorption spectroscopy[END_REF][START_REF] Kluczynski | Wavelength modulation absorption spectrometry -an extensive scrutiny of the generation of signals[END_REF] , offering a substantial decrease of 1/f noise, is a preferred technique in TDLAS to improve the sensitivity and robustness via transforming the detection frequency to a higher region. In contrast to direct absorption spectroscopy, the absolute concentration of a targeted gas cannot be determined from WMS spectral signals: the intensity of the WMS signal is dependent on pressure, temperature and laser output power etc. In order to lessen the common-mode noise, Rieker et al. proposed the 2f/1f normalization [START_REF] Rieker | Calibration-free wavelength-modulation spectroscopy for measurements of gas temperature and concentration in harsh environments[END_REF] , and Wang et al. developed a 2f/DS-sine to further simplify the calibration approach [START_REF] Wang | Laser frequency locking and intensity normalization in wavelength modulation spectroscopy for sensitive gas sensing[END_REF] . Moreover, Qu et al. presented a calibration-free scanned wavelength modulation spectrometer to derive species concentration and temperature by fitting to the 2f/1f-WMS lineshapes [START_REF] Qu | Calibration-free scanned wavelength modulation spectroscopy--application to H(2)O and temperature sensing in flames[END_REF] . The previous calibration-free WMS methods, in most cases, were based on a digital lock-in technique that demanded huge amounts of computing resources and was usually time-consuming. For the goal of fast response measurements, analog lock-in technology offers the function of instant demodulating process without occupying computing resources. Therefore, it is suitable for real-time WMS detection, especially for some cases including complex spectral line fitting procedures.

CH4 is of significant importance due to its strong greenhouse effect, and it's use in industry and urban utility, so its detection is critical for global warming and safety in human activity [START_REF] Willis | GHG Methodologies for Sewer CH4, Methanol-Use CO2, and Biogas-Combustion CH4 and their Significance for Centralized Wastewater Treatment[END_REF][START_REF] Saunois | The global methane budget 2000-2012[END_REF][START_REF] Vojtisek-Lom | On-road and laboratory emissions of NO, NO2, NH3, N2O and CH4 from late-model EU light utility vehicles: Comparison of diesel and CNG[END_REF] . On-board compact laser sensors for patrol vehicles can enormously increase the efficiency of leakage inspection. The strong fingerprint fundamental bands for CH4 locate at around 3.3 μm and 7.3 μm, so it is easy to reach sensitive detection in the corresponding regions. Based on the WMS, Qu et al. used mid-IR quantum cascaded laser (QCL) to detect CH4 at 7.5 μm and achieved a sensitivity of 0.04 ppm [START_REF] Qu | Mid-Infrared Trace CH4 Detector Based on TDLAS-WMS[END_REF] . Dong et al. monitored CH4 by using a 3.3 μm interband cascade laser and achieved minimum detection limit of 1.4 ppb for a 60 s averaging time [START_REF] Dong | Compact CH4 sensor system based on a continuouswave, low power consumption, room temperature interband cascade laser[END_REF] . Nevertheless, mid-IR diode lasers and detectors are more expensive, bulky and less durable compared to robust near-IR diode lasers. Thus, near infrared diode lasers are more appropriate for integration into small-sized laser sensors. Recently, Liu et al. [START_REF] Liu | Highly sensitive detection of methane by near-infrared laser absorption spectroscopy using a compact dense-pattern multipass cell[END_REF] and Dong et al. [START_REF] Cui | Three-Dimensional Printed Miniature Fiber-Coupled Multipass Cells with Dense Spot Patterns for ppb-Level Methane Detection Using a Near-IR Diode Laser[END_REF] , respectively, developed the novel multi-pass cells with low-volume and achieved a minimum detection limit of tens to hundreds ppb in the near-IR region.

In this work, a palm-sized spectrometer was developed based on TDLAS combining with a novel double-layer toroidal absorption cell (DLTAC). A home-made electronic circuit system was designed for laser control, signal acquisition, demodulation and data communication. Calibration-free scanned (CFS)-WMS was employed and a lineshape fitting procedure was carried out for deriving the concentration of the trace gas. The measurement of CH4 was performed by a distributed feedback (DFB) diode laser. The characteristics, such as response time and robustness towards temperature variation and mechanical vibration were measured for characterizing the performance of the spectral sensor. Allan deviation was utilized to analyze the stability of the sensor. To validate the application of the sensor, a demonstration of atmospheric CH4 monitoring was shown by two days of continuous detection.

Principle of CFS-WMS

For WMS, both the laser intensity and laser wavelength (frequency) are modulated by a high frequency sinusoidal modulation (fm) superposed in a low frequency sawtooth scan (fs). For one sawtooth period, the expression of the wavelength ν(t) is described as [START_REF] Qu | Calibration-free scanned wavelength modulation spectroscopy--application to H(2)O and temperature sensing in flames[END_REF] :

ν(t) = ν s (t) + a(t)cos(2πf m t + θ) , (1)
where νs(t) is the center laser wavelength without modulation and a(t) is the modulation depth that relies on the amplitude of the modulation signal, θ is the initial phase of wavelength modulation. According to Lambert-Beer's law, the laser intensity transmitted through the absorbing species is represented as:

I(t) = I 0 (t)exp[-α(ν)] = I 0 (t)exp[-cPLS(T)φ(ν)],
(2) where I0(t) is the laser intensity without absorbing medium and α(ν) is the absorbance defined by the concentration of the detected gas c, total pressure P, optical length L. Herein, line strength of transition S(T) relies on the temperature of the gas T in the cell and on the line shape function φ(ν). The harmonic signals Snf are extracted by a lock-in amplifier:

S nf = I(t)sin(2πnf m t) ⊗ LPF 2 + I(t)cos(2πnf m t) ⊗ LPF 2 , ( 3
)
where LPF is a low-pass filter. The WMS-2f signal is the preferred option relating with the signal to noise (SNR) and the signal amplitude at the center of the absorption line. The 2nd-harmonic signal normalized by the 1stharmonic signal (WMS-2f/1f) [START_REF] Salati | Apodized 2f/1f wavelength modulation spectroscopy method for calibration-free trace detection of carbon monoxide in the near-infrared region: theory and experiment[END_REF] can improve the robustness and reliability of sensors in harsh environments. The concentration of the absorbing gas is determined by least-square fitting of the observed WMS-2f/1f signal to lineshape function. For reliable fitting, temperature and pressure were measured by built-in sensors. The amplitude of WMS-2f/1f signal depends dramatically on both modulation amplitude (MA) and phase shift (PS) between laser intensity and wavelength modulation, therefore these factors need to be known before fitting WMS-2f/1f signal. In the present CFS-WMS scheme, MA and PS are determined by least-square fitting the direct absorption spectral signal to the measured signal in the presence of the absorbing gas, as depicted in Fig. 1(a). The temperature and pressure parameters in Eq.( 2) can be obtained by built-in sensor. Then, the concentration was derived by fitting the observed WMS-2f/1f signal using a Lorentzian function [START_REF] Topfer | Room-temperature mid-infrared laser sensor for trace gas detection[END_REF] , plotted in Fig. 1(b).

H2O and CO2 are two perturbing species to the CH4 spectral measurement near 1.65 µm, so it is necessary to choose a discrete and strong CH4 absorption line. Based on the HITRAN 2016 database [START_REF] Gordon | The HITRAN2016 molecular spectroscopic database[END_REF] , the absorption lines for 30 ppm CH4, 1% H2O and 0.3% CO2 around 6047 cm -1 are simulated at 1 atm, 300 K and 8.4 m absorption path, and are plotted in Fig. 2. As depicted in Fig. 2, the CH4 line is discret and the influence from H2O can be neglected. To ensure that the H2O concentration was less than 0.1% even for applications at high humidity, a drier was applied in our measurements.

Experimental

The schematic diagram of the palm-sized spectrometer for CH4 is shown in Fig. 3(a). The output of a fibercoupled distributed feedback (DFB) diode laser (NLK1U5FAAA, NEL) operated at 1.653 μm and room temperature passed through an optical collimator and was then coupled to a novel DLTAC. In principle, it is more suitable for ICL or QCL lasers emitting at the mid-infrared spectral region than the near infrared one due to our gold-coated inner surface with a higher reflectivity. The laser scanned from 6043 to 6050 cm -1 , covering the CH4 line at 6046.96 cm -1 . The temperature of the TEC was fine-tuned around 25.1 ℃ for guaranteeing the absorption peak of CH4 being located at the center of the scanning range. The output laser intensity from the DLTAC was collected by a homemade extended InGaAs detector (Netopto, NT52-1K-N). The thermometer and barometer sensors were used to monitor the parameters of the ambient atmosphere. As shown in Fig. 3 (b), the physical dimension of the sensor is minimized into 24×15×16 cm 3 due to the compact optical and home-made electronic design. The sensor is made of aluminum alloy, a low dense material, while the thin toroidal multipass cell is made of stainless steel. The total weight of the sensor is 2.2 kg. The wavelength modulation and WMS-2f/1f signal were applied to enhance the detection sensitivity.

The effective optical path length of the DLTAC is 8.35 m with 84 reflections exhibiting a polygon star spot pattern of two layers. It was indicated that the optical interference was weak, and the signal amplitude was rational for our designed toroid cell at 8 meters optical path length. The laser entrance and exit angles are both 8.57°. The number of reflections is determined by both the angle of incidence and the curvature radius in sagittal plane (Rs) and tangential plane (Rt). Herein, Rt is 99.57 mm and Rs is 50 mm. For illustrating the practical beam pattern, a 650 nm laser was used to form the visible beam trace, as presented in the insert of Fig. 3 (b).

The electronic part includes the temperature and current controller of the laser, the signal acquisition, temperature and pressure sensors, analog lock-in amplifier and data communication. The spectrometer is powered by business electricity with AC 220 V or by DC with 5 V for laser controller and 15 V for analog lock-in amplifier. The total electrical power consumption of the sensor is less than 10 W. The current sensor is powered by the homemade power board supplied by AC 220V, and it can also be powered by a battery. The laser controller module was designed on an 8 cm diameter PCB board. For the analog lock-in amplifier, the modulated signal at 7.572 kHz was generated by a voltage controlled oscillator (74HC4046A). A complex programmable logic device (CPLD, Altera, EPM7064AE) produced 1f and frequency-doubling (2f) signals and phase-shift. One of the 1f TTL signals from the CPLD output was smoothed by a 5th order Butterworth low-pass filter and superposed with slow sawtooth signal at 5 Hz. The output of the photo detector was first acquired by the ADC for determining the parameters of the laser and then demodulated by two AD630 lock-in amplifiers for 1f and 2f. Two LPFs were used for further filtration, respectively. The demodulated signals were collected by the ADC (AD7616) and then transmitted to the microcontroller unit (STM32F405). ADC was used for fitting the 2f/1f spectral signal transmitted by a serial port. The concentration of the trace gas can be derived by a least-square fitting procedure.

Results and discussion

In wavelength modulation spectroscopy, the amplitude of the WMS-2f/1f signal depends on the modulation amplitude [START_REF] Liu | Highly sensitive detection of methane by near-infrared laser absorption spectroscopy using a compact dense-pattern multipass cell[END_REF][START_REF] Sur | Multi-species laser absorption sensors for in situ monitoring of syngas composition[END_REF] . To determine the optimum modulation amplitude, WMS-2f signals of 50 ppm CH4 were measured at a series of modulation amplitudes at atmospheric pressure. From the measured signal, the maximum signal for CH4 is obtained at 0.27 V and the corresponding modulation index is around 2.2 [START_REF] Sur | Multi-species laser absorption sensors for in situ monitoring of syngas composition[END_REF][START_REF] Ste´ Phane Schilt | Wavelength modulation spectroscopy: combined frequency and intensity laser modulation[END_REF] , which is in accordance with the theoretical prediction. Thus, a modulation amplitude of 0.27 V was employed in the following measurements. To verify the CFS-WMS principle, the measured and simulated WMS-2f/1f signals of a series of different CH4 concentrations are presented in Tests were performed to investigate the performance of our sensor to mechanical vibrations and changes in ambient temperature. The vibration test is a crucial index to evaluate an on-board or airborne sensor. The toroidal multipass cell is constructed by a monolithic ring, so its instinct is immune to vibration effect. The slight intensity fluctuation of laser is mainly arisen from a tiny displacement between the incident beam and the toroidal multipass cell due to mechanical vibration. As presented in Figure 5(a), the vibration test was simulated by a shaker table (Fang Ce Precision instrument technology co. LTD, LD-35BL) with frequency from 0 to 400 Hz. As shown in Fig. 5(b), the concentration of the cell was set to 20 ppm. Overall, the difference between the measured and the real concentration was small in the entire range from 0 to 400 Hz. The largest error bar occurs at 200 Hz. It can be deduced that the natural frequency of our spectrometer is approximately 200 Hz. The test shows that the standard deviation of measurement is 0.36 ppm. The spectral intensity of peak-peak is sensitive to temperature variation, so one needs to either keep the temperature constant or correct the temperature influence if the sample concentration is determined by the spectral peak-peak amplitude. As depicted in Fig. 5(c), a test of our sensor to variation of the temperature is carried out by a high-low temperature test box (Wuxi yi'er da test equipment manufacturing Co., Ltd, GDW-150). The temperature interval was set to 10 ℃ and each temperature took 30 min to ensure the stability and reliability of the measurement. It illustrates that the peak-peak values of the spectra show an approximately linear decrease as the temperature increases from -10 to 50 ℃, as presented in Fig. 5(d). However, the temperature influence is diminished by line-fitting considering the temperature factor. The standard error, as result of temperature variation, is approximatly 0.11 ppm from -10 to 50 ℃. To evaluate the detection precision and the stability of the sensor, the signal of 20 ppm CH4 sealed in the cell was measured for about one hour. 34320 data points were obtained with the scan frequency of 5 Hz. As shown in Fig. 6(a), the Allan variance is a function of the measurement time. The figure indicates that the minimum detection limit is reached at 57.8 s and the minimum detectable concentration can be down to 22 ppb. The histogram of the measured concentration of CH4, helpful to evaluate the measurement precision, is plotted in Fig. 6(b). The data distribution resembles a Gaussian profile, red line in Fig. 6(b). The measurement precision, determined by the half width at half maximum (HWHM), was found to be 0.8 ppm and the R 2 value and the standard deviation value (σ) were calculated to be 0.986 and 0.680 ppm, respectively. The high sensitivity and high measurement precision confirm the high performance of the developed CH4 sensor based on toroidal multipass cell. The dynamic response of the sensor has been performed at the flow rate of 500 mL/min while the volume of the cell is about 157 mL. In Fig. 7(a), the rise time (10-90 %) and the fall time (90-10 %) are 16 s and 14 s, respectively. Fig. 7 By real-time and in-situ measurement of CH4 in ambient air for two days, the performance of the developed CH4 sensor was validated. The test was conducted on the campus of the Taiyuan University of Science and Technology. It is a mixed-use area with woods, apartments and nearby the Xizhonghuan Road with heavy traffic especially during rush hours. The temporal resolution of the sensor is determined by the scanning frequency of the diode laser (2 s) and the interval of gas extracted into the multi-pass cell (~20 s) which is usually shorter than the time to reach equilibrium of gases. Therefore, the time interval between the measurement and the process of extracting gas into the multi-pass cell does not affect the real-time measurement. During the measurement, a dryer and a particle filter were placed before the inlet of the multi-pass cell to remove humidity and aerosols. The sensor operated continuously during the measurement. The concentrations of ambient CH4 are shown in Fig. 8, while the data were captured at 0.5 Hz and averaged over 8 laser scans. The results of consecutive measurements in two days indicate the same CH4 concentration variation trends. The highest concentration of CH4 occurs at ~2am and the lowest concentration appears after 5am on both days. Moreover, the CH4 concentration is higher during the night than at daytime. The photochemical removal of CH4 from the atmosphere is very slow, leading to an atmospheric lifetime of CH4 of more than 10 years. Therefore, the variation in the CH4 concentrations during the day are possibly due to increased mixing of local CH4 sources after sunrise (6:40 am) while the increase in concentration after sunset (5:40 pm) is due to a decreased boundary layer height.

Conclusions

A portable optical sensor for trace gas quantification with a toroidal cell based on near infrared TDLAS is reported. The 8.35 m optical path length is realized with 84 reflections in two layers in a toroidal cell with an aspheric surface. The WMS-2f/1f calibration free method is employed to derive the concentration of the absorbing gas. The laser driver is implemented in a home-made electronic system. The minimum achieved detection limit of CH4 was 22 ppb with an integration time of 57.8 s. To validate the performance of the sensor, we continuously monitored the concentration of CH4 in ambient air for two days using a DFB laser emitting at 1.653 μm. The results indicate that the sensor is appropriate for monitoring trace gas.

Fig. 1 (Fig. 2 .

 12 Fig.1 (a) Measured (black dash line) and fitted simulated (red solid line) transmitted laser intensity; (b) Measured (black dash line) and fitted simulated (red solid line) WMS-2f/1f signal of 40 ppm CH4.

Fig. 3 .

 3 Fig. 3. (a) Schematic diagram of the experimental setup for CH4 detecting; (b) photograph of the real instrument

Figure 4 :Fig. 4 .

 44 Fig. 4. Linear relationship between the real concentrations and the fitted ones.

Fig. 5

 5 Fig.5 (a) Layout of the vibration test; (b) measured results obtained for 20 ppm CH4 at different vibration frequencies; (c) structure of the temperature test; (d) observed amplitude of peak-peak of WMS-2f/1f and best fitted concentration of 20 ppm CH4 at different temperatures.

Fig. 6 .

 6 Fig.6. (a) Row data of individual concentration measurements over 400 s (top) and Allan variance as the function of integration time

Fig. 7 (Fig. 8 .

 78 Fig.7 (a) Response time by varying the CH4 concentration between 0 and 20 ppm; (b) measured WMS-2f /1f CH4 signals at different concentrations.

Table of Content

 of Graphic For Table of Contents Only

Acknowledgement

The work presented in this paper was supported by National Natural Science Foundation of China (U1810129, 52076145, 11904252 & 42077201), State Key Laboratory of Applied Optics (SKLAO-201902), Transformation of Scientific and Technological Achievements Fund of Shanxi Province (201904D131025), Excellent Youth Academic Leader in Higher Education of Shanxi Province (2018), Key Research and Development Program of Shanxi Province of China (201803D31077, 201803D121090), Shanxi "1331 Project" Key Innovative Research Team (1331KIRT), Natural Science Foundation of Shanxi Province (201801D221017) and Fund for Shanxi Key Subjects Construction.


Philippe Castagliola
email: philippe.castagliola@univ-nantes.fr

Giovanni Celano
email: giovanni.celano@unict.it

Dorra Rahali
email: dorra.rahali@centralelille.fr

Shu Wu

Control Charts for Monitoring Time-Between-Events-and-Amplitude Data

Introduction

Today, due to the large availability of data, various kinds of processes can (and have to) be monitored using Statistical Process Monitoring (SPM) techniques based on advanced control charts. These kinds of processes can be of course industrial ones but, they can also be non industrial ones like in the biological/health-care (diseases, like the Covid-19 for instance), the geological (earthquakes or volcanic eruptions) or the accidental (traffic accidents, forest fires) fields. In all of these situations, people are usually focusing on two characteristics:

1. the time T between two consecutive specific (usually, adverse) events of interest E, 2. the amplitude X of each of these events. The characteristics T and X defined above are the key factors to be monitored for an event and they are usually referred to as the TBEA (Time Between Events and Amplitude) characteristics. In general, a decrease in T and/or an increase in X can result in a negative, hazardous or disastrous situation that needs to be efficiently monitored with control charts.

The first TBE (i.e. without taking into account the amplitude characteristic) type of control chart goes back to [START_REF] Calvin | Quality Control Techniques for "Zero Defects[END_REF], who proposed to monitor the cumulative number of conforming items between two non-conforming ones. The initial idea was to find a method to improve the traditional attribute control charts that are known to be ineffective in the case of high-quality processes in which the occurence of non-conforming products is very rare. This initial idea has then been investigated by [START_REF] Lucas | Counted Data CUSUM's[END_REF] and [START_REF] Vardeman | Average Run Lengths for CUSUM Schemes when Observations are Exponentially Distributed[END_REF] and, subsequently, many other researchers started to contribute to this area. [START_REF] Radaelli | Planning Time-Between-Events Shewhart Control Charts[END_REF] proposed to design and implement oneand two-sided Shewhart-type TBE control charts assuming that the counts can be modeled as a homogeneous Poisson process. [START_REF] Gan | Designs of One-and Two-Sided Exponential EWMA Charts[END_REF] developed an EWMA (Exponentially Weighted Moving Average) control chart monitoring the rate of occurrences of rare events based on the inter-arrival times of these events. [START_REF] Benneyan | Number-Between g-Type Statistical Quality Control Charts for Monitoring Adverse Events[END_REF] used the geometric ("g" chart) and the negative binomial (called "h" chart) distributions in order to monitor the number of cases between hospital-acquired infections. [START_REF] Xie | Some Effective Control Chart Procedures for Reliability Monitoring[END_REF] proposed a control chart for TBE data based on the exponential distribution while [START_REF] Borror | Robustness of the Time Between Events CUSUM[END_REF] extended it using a CUSUM (Cumulative Sum) scheme and evaluated its robustness in the case of the Weibull and lognormal distributions. [START_REF] Liu | A Comparative Study of Exponential Time Between Events Charts[END_REF] compared the ATS (Average Time to Signal) performance of several continuous TBE charts including the CQC, CQC-r, exponential EWMA and exponential CUSUM charts. Zhang et al. (2007) investigated the case of gamma distributed TBE data and they developed a control chart based on a random-shift model to compute the out-of-control ATS. In the case of multistage manufacturing processes, [START_REF] Shamsuzzaman | Integrated Control Chart System for Time-Between-Events Monitoring in a Multistage Manufacturing System[END_REF] developed a control chart for TBE data and designed it using a statistical oriented approach while Zhang et al. (2011a) designed it using a first economic oriented approach and Zhang et al. (2011b) developed it using a second economic oriented approach assuming random process shifts. The use of supplementary runs rules has also been proposed for monitoring TBE data by [START_REF] Cheng | An ARL-Unbiased Design of Time-Between-Events Control Charts with Runs Rules[END_REF]. [START_REF] Qu | Time-Between-Event Control Charts for Sampling Inspection[END_REF] studied some TBE control charts that can be used for sampling inspection. [START_REF] Shafae | Cumulative Sum Control Charts for Monitoring Weibull-Distributed Time Between Events[END_REF] evaluated the performance of three TBE CUSUM charts and [START_REF] Fang | Monitoring of Time-Between-Events with a Generalized Group Runs Control Chart[END_REF] proposed a generalized group runs TBE chart for a homogenous Poisson failure process.

The first paper that proposed a combined scheme for monitoring the time interval T of an event E as well as its amplitude X has been introduced by [START_REF] Wu | A Control Scheme for Monitoring the Frequency and Magnitude of an Event[END_REF] who refered it to as a TBEA (Time Between Events and Amplitude) chart. After this paper, several single TBEA charts have been developed, see for instance [START_REF] Qu | A CUSUM Scheme for Event Monitoring[END_REF], [START_REF] Cheng | Simultaneously Monitoring Frequency and Magnitude of Events Based on Bivariate Gamma Distribution[END_REF], [START_REF] Ali | Time and Magnitude Monitoring based on the Renewal Reward Process[END_REF], [START_REF] Qu | A CUSUM Chart for Detecting the Intensity Ratio of Negative Events[END_REF] and, very recently, [START_REF] Sanusi | Simultaneous Monitoring of Magnitude and Time-Between-Events Data with a Max-EWMA Control Chart[END_REF].

As it can be noticed, this stream of research is rather recent and few publications have already been devoted to. Therefore, the goal of this chapter is to further investigate it and to hopefully open new research directions. More specifically, this chapter will be splitted into three parts:

1. In section 2 we will introduce and compare three different statistics, denoted as Z 1 , Z 2 and Z 3 , suitable for monitoring TBEA data, in the case of four distributions (gamma, lognormal, normal and Weibull), when the time T and the amplitude X are considered as independent random variables.

2. In section 3, we will compare the three statistics introduced in section 2, for the same distributions, but considering that the time T and the amplitude X are dependent random variables. A model based on three types of Copulas will be used to define the dependence between T and X. 3. Finally, in section 4, in order to overcome the "distribution choice" dilemma, we will introduce a distribution-free approach coupled with an upper-sided EWMA scheme. In addition, a specific technique called "continuousify" will be presented in order to compute the Run Length properties of the proposed upper-sided EWMA TBEA control chart in a reliable way.

2 TBEA charts for independent Times and Amplitudes

Model

Let D 0 = 0, D 1 , D 2 , . . . be the dates of occurrence of a specific negative event E, let T 1 = D 1 -D 0 , T 2 = D 2 -D 1 , . . . be the time intervals between two consecutive occurrences of the event E and let X 1 , X 2 , . . . be the corresponding magnitudes of this event occurring at times D 1 , D 2 , . . . (see Figure 1). It must be noted that D 0 = 0 is the date of a "virtual" event which has no amplitude associated with. In this section, we assume that T and X are two mutually independent continuous random variables, both defined on [0, +∞). Let F T (t|θ T ) and F X (x|θ X ) be the c.d.f. (cumulative distribution function) of T and X, respectively, and let f T (t|θ T ) and f X (x|θ X ) be the p.d.f. (probability distribution function) of T and X, respectively, where θ T and θ X are the corresponding vector of parameters. Let also define µ T = E(T ), µ X = E(X), σ T = σ(T ) and σ X = σ(X) be the expectation and standard-deviation of T and X, respectively. By definition, when the process is in-control, we have

D0 = 0 D1 D2 D3 • • • • • • • • • T1 T2 T3 X1 X2 X3 Times      Dates      Amplitudes             
θ T = θ T0 , θ X = θ X0 , µ T = µ T0 , µ X = µ X0 , σ T = σ T0 , σ X = σ X0 and, when the process is out-of-control, we have θ T = θ T1 , θ X = θ X1 , µ T = µ T1 , µ X = µ X1 , σ T = σ T1 , σ X = σ X1 .
Because the reference scales for the random variables T and X can be very different and, in order to not favour one random variable over the other one, we suggest to define (and work with) the "normalized to the mean" new random variables T and X as the in-control standardized counterparts of T and X, i.e.

T = T µ T0 , X = X µ X0 .
Clearly, when the process is in-control we have E(T ) = E(X ) = 1.

Statistics to be monitored

In order to simultaneously monitor the time T between an event E and its amplitude X, we suggest to define several dedicated statistics Z = Z(T , X ), functions of the random variables T and X , satisfying the following two properties:

Z ↑ if either T ↓ or X ↑, (1)
Z ↓ if either T ↑ or X ↓ . (2)
Of course, there are many possible choices for the statistic Z. A first possible choice for the statistic Z (denoted as the Z 1 statistic) satisfying properties (1) and ( 2) is simply

Z 1 = X -T . (3)
This random variable is defined on (-∞, +∞) and its c.d.f. F Z1 (z|θ Z ) and p.d.f. f Z1 (z|θ Z ) can be obtained by integrating (see Figure 2 (a) and (b)) over all the couples (X , T ) ∈ R +2 satisfying Z 1 = X -T ≤ z, and they are equal to where θ Z = (θ T , θ X ) is the corresponding combined vector of parameters.

F Z1 (z|θ Z ) = 1 -µ X0 +∞ 0 F T ((x -z)µ T0 |θ T )f X (xµ X0 |θ X )dx, (4)
f Z1 (z|θ Z ) = µ T0 µ X0 +∞ 0 f T ((x -z)µ T0 |θ T )f X (xµ X0 |θ X )dx, ( 5
) X -T ≤ z X T 0 X -T ≤ z X T 0 z (a) (b) X T ≤ z X T 0 X + 1 T ≤ z z X T 0 (c) (d)
A second possible choice for the statistic Z (denoted as the Z 2 statistic) satisfying properties (1) and ( 2) is defining it as the ratio between the two characteristics of an event E:

Z 2 = X T . (6)
This random variable is defined on [0, +∞) and its c.d.f. F Z2 (z|θ Z ) and p.d.f. f Z2 (z|θ Z ) can be obtained by integrating (see Figure 2 (c)) over all the couples (X , T ) ∈ R +2 satisfying Z 2 = X T ≤ z, and they are equal to

F Z2 (z|θ Z ) = 1 -µ X0 +∞ 0 F T xµ T0 z |θ T f X (xµ X0 |θ X )dx, (7) f Z2 (z|θ Z ) = µ T0 µ X0 z 2 +∞ 0 xf T xµ T0 z |θ T f X (xµ X0 |θ X )dx. (8)
Finally, a third possible choice for the statistic Z (denoted as the Z 2 statistic) satisfying properties (1) and ( 2) is

Z 3 = X + 1 T . (9)
This random variable which should be considered as a hybrid of the two previous ones is also defined on [0, +∞) and its c.d.f. F Z3 (z|θ Z ) and p.d.f. f Z3 (z|θ Z ) can be obtained by integrating (see Figure 2 (d)) over all the couples (X , T ) ∈ R +2 satisfying Z 3 = X + 1 T ≤ z, and they are equal to

F Z3 (z|θ Z ) = F X (zµ X0 |θ X ) -µ X0 z 0 F T µ T0 z -x |θ T f X (xµ X0 |θ X )dx, (10)
f Z3 (z|θ Z ) = µ T0 µ X0 z 0 1 (z -x) 2 f T µ T0 z -x |θ T f X (xµ X0 |θ X )dx. ( 11
)
More details on how to derive the c.d.f. and p.d.f. of statistics Z 1 , Z 2 and Z 3 provided above can be found in the Appendix section of [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Several Distributions[END_REF]. Concerning these c.d.f. and p.d.f. it has to be noted that it is generally not possible to obtain a closed form solution for them and the only solution is to numerically compute these ones by using quadrature techniques.

Control limit

As it is more important to detect an increase in Z (in order to avoid more damages or injuries, for instance) rather than a decrease, we suggest to only define an upper control limit UCL Z for the TBEA charts based on statistics

Z ∈ {Z 1 , Z 2 , Z 3 } as UCL Z = F -1 Z (1 -α|θ Z0 ), ( 12
)
where α is the type I error, θ Z0 = (θ T0 , θ X0 ) and F -1 Z (. . . |θ Z0 ) is the inverse c.d.f. of Z numerically obtained by solving equation F Z (z|θ Z0 ) = α for z using a one dimension root finder.

Time to Signal properties

The type II error β of the upper-sided TBEA charts based on statistic

Z ∈ {Z 1 , Z 2 , Z 3 } is equal to β = F Z (UCL Z |θ Z1 ), (13)
where θ Z1 = (θ T1 , θ X1 ). The out-of-control ATS (Average Time to Signal) and SDTS (Standard Deviation Time to Signal) of the upper-sided TBEA charts based on statistic Z ∈ {Z 1 , Z 2 , Z 3 } can be obtained using the expectation and variance of compound random variables (see also [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Several Distributions[END_REF] for more details) and they are equal to

ATS 1 = µ T1 1 -β , ( 14
) SDTS 1 = σ 2 T1 1 -β + µ 2 T1 β (1 -β) 2 . ( 15
)
When the process is in-control, we have 1β = α and, consequently, we have the following equivalence for the in-control ATS

ATS 0 = µ T0 α ⇔ α = µ T0 ATS 0 .

Comparative studies

As in [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Several Distributions[END_REF], in order to compare the three TBEA charts defined in sub-section 2.2 and based on the statistics Z ∈ {Z 1 , Z 2 , Z 3 } we have chosen to investigate four different types of distribution that are only dependent on two parameters a and b. The choice of the Gamma, Lognormal, Normal and Weibull distributions is driven by the fact that these ones are very often selected to model time oriented random variables. For this reason, the two parameters beta distribution has been excluded from the benchmark as it is rarely selected for representing time oriented variables. Of course, more complex distributions could have been considered (like the four parameters Beta or Johnson's distributions) but, due to the fact that only the nominal mean µ 0 and standard-deviation σ 0 are assumed to be known, we restricted our choice to a selection of two parameters distributions. These distributions are summarized in Table 1 with their names, parameter settings and p.d.f. f (x|a, b).

In this table, f Nor (. . . ) stands for the p.d.f. of the normal (0, 1) distribution.

Gamma a > 0, b > 0 exp(-x b )x a-1 b a Γ (a) Lognormal a, b > 0 b x fNor(a + b ln(x)) Normal a, b > 0 1 b fNor x-a b Weibull a > 0, b > 0 a b x b a-1 exp -x b a
More specifically, we have selected 9 different in-control configurations to be investigated for T (the normal distribution has not been considered as a possible choice for the time between events) and 11 in-control configurations to be investigated for X, i.e. a total of 99 scenarios for (T, X). All of them (see Table 2) are such that the in-control mean is µ 0 = 10 and the in-control standard-deviation is σ 0 ∈ {1, 2, 5} (except for the normal distribution, where σ 0 ∈ {1, 2} only). In addition, Table 2 also provides the values of the in-control parameters a 0 and b 0 and the corresponding skewness coefficient γ 0 . The upper control limits UCL Z for the three TBEA charts based on statistics Z ∈ {Z 1 , Z 2 , Z 3 }, satisfying ATS 0 = 370.4, have been obtained in Table 3 for the 99 possible scenarios defined in Table 2. As it can be seen, regardless of the statistic Z and the distribution of X, these upper control limits tend to be similar if σ 0 is small (say σ 0 = 1) and the distribution of T is either gamma or lognormal. But, when T follows a Weibull distribution, the control limits are larger than those of the gamma or lognormal distributions.

When an out-of-control situation occurs in a TBEA process (corresponding to an upper shift in Z), it can be due to i) a mean shift only in the time T from µ T0 to µ T1 = δ T µ T0 , ii) a mean shift only in the amplitude X from µ X0 to µ X1 = δ X µ X0 , or iii) a combination of the two previous cases, where δ T ≤ 1 and δ X ≥ 1 are the parameters quantifying the change in the time and amplitude, respectively. But, as the actual values of δ T and δ X are usually unknown by the practitioner, it is therefore difficult to evaluate the three TBEA charts based on statistics Z ∈ {Z 1 , Z 2 , Z 3 } using the ATS 1 criterion defined in ( 14) that depends on a specific values for δ T and/or δ X . For this reason, it is therefore preferable to use the following more general criterion denoted as EATS 1 (Expected Average Time to Signal) and defined as

EATS 1 = δ T ∈Ω T δ X ∈Ω X f δ T (δ T )f δ X (δ X )ATS 1 (δ T , δ X ),
where Ω T and Ω X are the sets of the potential shifts for δ T and δ X , respectively, and f δ X (δ X ) and f δ T (δ T ) are the probability mass functions of the shifts δ T and δ X over the sets Ω T and Ω X . In this chapter, we adopt the classical assumption that confines f δ T (δ T ) and f δ X (δ X ) to be discrete uniform distributions over Ω T and Ω X , respectively. If we want to investigate i) a mean shift only due to the time T then we suggest to fix Ω T = {0.5, 0.55, . . . , 0.9, 0.95} and Ω X = {1}, ii) a mean shift only due to the amplitude X then we suggest to fix Ω T = {1} and Ω X = {1.1, 1.2, . . . , 1.9, 2} and iii) a mean shift due to the time T and the amplitude X then we suggest to fix Ω T = {0.5, 0.55, . . . , 0.9, 0.95} and Ω X = {1.1, 1.2, . . . , 1.9, 2}.

For the 99 possible scenarios defined in Table 2, Table 4 gives the EATS 1 values of the three TBEA charts based on statistics Z ∈ {Z 1 , Z 2 , Z 3 } when Ω T = {0.5, 0.55, . . . , 0.9, 0.95} and Ω X = {1.1, 1.2, . . . , 1.9, 2}. Values of EATS 1 in bold characters are the smallest ones among statistics Z 1 , Z 2 or Z 3 . From Table 4, it can be deduced that when there is a shift in both T and X, the most efficient statistic is Z 1 (in 56% of the cases with an average EATS 1 value EATS 1 = 14.91), followed by statistic Z 2 (in 35% of the cases with EATS 1 = 29.15) and, finally, statistic Z 3 (in only 5% of the cases with EATS 1 = 11.74).

The cases where the mean shift is only due to the time T or the mean shift is only due to the amplitude X have both been investigated in Rahali et al.

(2019) (see Tables 3 and4, pages 245-246). In these cases, the conclusions are

• if the mean shift is only due to the time T , then the most efficient statistic is Z 3 (in 71% of the cases with EATS 1 = 59.88) followed by Z 2 (in 29% of the cases with EATS 1 = 71.82) while the statistic Z 1 never provides the smallest EATS 1 and should not be considered here as an efficient monitoring statistic. • if the mean shift is only due to the amplitude X, then the statistic Z 1 is the best option as it always gives the smallest EATS 1 values, regardless of the combination under consideration. In this case, Z 2 and Z 3 should not be considered as potential efficient monitoring statistics.

Illustrative example

This illustrative example has been detailed for the first time in [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Several Distributions[END_REF] and it is based on a real data set concerning the time (T in days)

Table 4. EATS1 values when ΩT = {0.5, 0.55, . . . , 0.9, 0.95} and ΩX = {1.1, 1.2, . . . , 1.9, 2} for the three TBEA charts based on statistics Z ∈ {Z1, Z2, Z3} 19.0 20.3 38.3 19.0 20.6 42.4 19.0 19.9 18.9 19.7 33.7 73.3 73.3 73.5 73.3 73.3 73.9 73.3 73.3 73.3 73.3 73.5 between fires in forests of the region "Provence -Alpes -Côte D'Azur" in the south-east of France and their amplitudes (X measured as the burned surface in ha = 10000m 2 ). This data set reports a total of 92 significant fires from 2016/10 to 2017/9: the data set has been split into m = 47 fires occuring during the "low season" from 2016/10 to mid 2017/6 (used as Phase 1 data) and n = 45 fires occuring during the "high season" from mid 2017/6 to 2017/9 (used as Phase 2 data). The values of T and X are recorded in Table 5 (as well as the values of the statistics Z 1 , Z 2 and Z 3 ) and they are also plotted in Figure 3 where it is clear that the values of T during the "high season" are smaller than those during the "low season" and the values of X during the "high season" are larger than those during the "low season".

Statistic Z1 T X → Gamma Lognormal Normal Weibull ↓ σ0 1 2 5 1 2 5
The use of the Kendall's and Spearman's rank correlation tests on the whole data set yields p-values larger than the significance level of 0.05 (0.2 for the Kendall's test and 0.19 for the Spearman's test) validating the fact that the random variables T and X are uncorrelated (a key assumption in this section). Among the four distributions considered in Table 1, the use of the Kolmogorov-Smirnov's test shows that the best fit for both T and X is the lognormal distribution with parameters (a 0 = -1.2648, b 0 = 1.0302) for T and (a 0 = -1.6697, b 0 = 0.8624) for X.

The three TBEA charts, corresponding to the statistics Z ∈ {Z 1 , Z 2 , Z 3 } are plotted in Figure 4 along with their upper control limits UCL Z1 = 6.0306, UCL Z2 = 28.1209 and UCL Z3 = 19.3885 (assuming ATS 0 = 730, i.e. 2 years).

As it can be seen, these charts detect several out-of-control situations during the "high season" confirming that a decrease in the time between fires occured as well as an increase in the amplitude of these fires.

3 TBEA charts for dependent Times and Amplitudes

Motivation

In the previous section, the time T between events and their amplitudes X have been considered as independent random variables. But, in practice, this is not always the case. For example, there are natural situations for which the amplitudes tend to become larger when the times between events become shorter (i.e. a negative correlation). Such kind of situations is likely to occur for example in the case of earthquakes for which, in a first phase, small amplitude earthquakes may occur with a low frequency (large time between events) and, suddenly, in a second phase, the occurence frequency of these earthquakes may increase (shorter time between events) with a negatively correlated increase in their amplitudes. The same kind of situations may also arise in the case of forest fires occuring, in a first phase during the "humid season", with a low frequency and small amplitudes (surfaces burned) and Sep 17

Phase 1 Phase 2 X Fig. 3. Time (T in days) between fires and amplitudes (X as the burned surface in ha) corresponding to the data set in Table 5 becoming more disastrous during the "dry season" with shorter time between the occurrence of forest fires and larger amplitudes (see for instance the 2019 forest fires in Amazonia or Siberia or the 2020 forest fires in Australia and USA). Positive correlation between T and X (i.e. the time between events becomes shorter and the amplitude becomes smaller) is also possible as the the forthcoming illustrative example in this section will depict it.

Until now, very few research papers have investigated TBEA control charts by considering the potential dependence between the two variables T and X.

- Cheng and Mukherjee (2014) were the first to investigate a T 2 TBEA control chart by using a bivariate SAT (Smith-Adelfang-Tubbs) Gamma distribution to model the joint probability of T and X. This work has been extended later by [START_REF] Cheng | Simultaneously Monitoring Frequency and Magnitude of Events Based on Bivariate Gamma Distribution[END_REF] who developed a similar approach based on a MEWMA (multivariate exponentially weighted moving average) procedure.

In this section, instead of specifying a particular bivariate joint distribution for (T, X), like the SAT distribution for instance, we will assume that the marginal distributions of T and X are both known (and they can be almost anything) and we will use the Copulas mechanism (popularized by [START_REF] Sklar | Fonctions de Répartition à n Dimensions et leurs Marges[END_REF]) in order to model the dependence between the time T and the amplitude X. The use of Copulas in the Statistical Process Monitoring field is not so common. We can cite for instance [START_REF] Fatahi | A Bivariate Control Chart Based on Copula Function[END_REF], [START_REF] Dokouhaki | A Copula Markov CUSUM Chart for Monitoring the Bivariate Auto-Correlated Binary Observations[END_REF], [START_REF] Busababodhin | Copula Modelling for Multivariate Statistical Process Control: A Review[END_REF] and Sukparungsee et al. ( 2018) who all proposed various types of control charts based on Copulas.

Model

In this section, we assume that (X, T ) ∈ R 2 + and their joint continuous c.d.f. is equal to

F (T,X) (t, x|θ T , θ X , θ) = C(F T (t|θ T ), F X (x|θ X )|θ), (16)
where F T (t|θ T ) and F X (x|θ X ), as defined in section 2.1, are the marginal c.d.f. of T and X, respectively, C(u, v|θ) is a Copula containing all information on the dependence structure between T and X, and θ is a dependence parameter quantifying the dependence between the marginals. In addition, let

f (T,X) (t, x|θ T , θ X , θ) = c(F T (t|θ T ), F X (x|θ X )|θ)f T (t|θ T )f X (x|θ X ) (17)
be the joint p.d.f. of (X, T ) where f T (t|θ T ) and f X (x|θ X ) are the marginal p.d.f.'s of T and X, respectively, and c(u, v|θ) = ∂C(u,v|θ) ∂u∂v is the Copula density. As explained in section 2.1, in order to not favor one random variable over the other one, the new random variables T = T µ T 0 and X = X µ X 0 are introduced as the in-control standardized counterparts of T and X, where µ T0 and µ X0 are the (known) in-control expectation of T and X, respectively. It is easy to show that the joint c.d.f. F (T ,X ) (t, x|θ T , θ X , θ) and joint p.d.f.

f (T ,X ) (t, x|θ T , θ X , θ) of (X , T ) ∈ R 2 + are equal to F (T ,X ) (t, x|θ T , θ X , θ) =C(F T (tµ T0 |θ T ), F X (xµ X0 |θ X )|θ), f (T ,X ) (t, x|θ T , θ X , θ) =µ T0 µ X0 c(F T (tµ T0 |θ T ), F X (xµ X0 |θ X )|θ) f T (tµ T0 |θ T )f X (xµ X0 |θ X ).
The closed form formulas for the c.d.f. and p.d.f. of the statistics Z 1 , Z 2 and Z 3 defined in section 2.2 are provided in [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Correlated Data[END_REF]. The definition of the upper control limit UCL Z of the TBEA charts with dependent times T and amplitudes X is similar to the one in equation ( 12) and it just requires the addition of the dependence parameter θ, i.e. UCL Z = F -1 Z (1 -α|θ Z0 , θ). The formulas for computing ATS 1 and SDTS 1 are the same as in equations ( 14) and (15), respectively.

Comparative studies

In order to compare the three TBEA charts for dependent times T and amplitudes X based on statistic Z ∈ {Z 1 , Z 2 , Z 3 }, the same distributions listed in Table 1 have been chosen and the same 99 possible scenarios defined in Table 2 have been investigated. The Archimedean bivariate Copulas of [START_REF] Gumbel | Distributions des Valeurs Extremes en Plusieurs Dimensions[END_REF], [START_REF] Clayton | A Model for Association in Bivariate Life Tables and its Application in Epidemiological Studies of Familial Tendency in Chronic Disease Incidence[END_REF] and [START_REF] Frank | On the Simultaneous Associativity of F (x, y) and x + y -F (x, y)[END_REF] have been chosen in this section to model the dependence between T and X. The Gumbel's (also called Gumbel-Hougaard's) and Clayton's Copulas are two asymmetric Copulas exhibiting a larger dependence in the positive tail than in the negative one (for the Gumbel's Copula) and in the negative tail than in the positive one (for the Clayton's Copulas). The Frank's Copula is a symmetric one that can be used to model dependence structures with either positive or negative correlation. Other Archimedean bivariate Copulas could have also been investigated, like for instance the [START_REF] Ali | A Class of Bivariate Distributions including the Bivariate Logistic[END_REF] and [START_REF] Joe | Parametric Families of Multivariate Distributions with Given Margins[END_REF] Copulas but, for simplicity and also due to their popularity, we only confined our investigations to the Gumbel's, Clayton's and Frank's Copulas. Details concerning the definition of each of these Copulas C(u, v|θ), the domain of definition for θ and the relationship between the Kendall's rank correlation coefficient τ and the dependence parameter θ are provided in Table 6. In order to facilitate the use of these Copulas, Table 7 simply provides pre-computed values of θ for several selected values of the Kendall's rank correlation coefficient τ ∈ {0, 0.1, 0.2, . . . , 0.9}.

When it is not possible to model a negative dependence with the Copulas defined above, it is possible to use 90 or 270 degrees rotated versions C 90 or C 270 of these Copulas using the following transformations (see [START_REF] Brechmann | Modeling Dependence with C-and D-Vine Copulas: The R Package CDVine[END_REF])

C 90 (u 1 , u 2 ) = u 2 -C(1 -u 1 , u 2 ), C 270 (u 1 , u 2 ) = u 1 -C(u 1 , 1 -u 2 ).
Similarly to Table 3 (for independent times T and amplitudes X) and assuming ATS 0 = 370.4, the upper control limits UCL Z of the three TBEA charts with dependent times T and amplitudes X are reported in Tables 345of [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Correlated Data[END_REF] for the 99 scenarios defined in Table 2 and for the 3 Copulas defined above. From Tables 3-5 in [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Correlated Data[END_REF], the following conclusions can be drawn

) θ + (-ln(v)) θ 1 θ [1, ∞) τ = 1 -1 θ ⇔ θ = 1 1-τ Clayton max(0, u -θ + v -θ -1) -1 θ [-1, ∞)\{0} τ = θ θ+2 ⇔ θ = 2τ 1-τ
Frank -1 θ ln 1 + (e -θu -1)(e -θv -1) • For a fixed statistic Z ∈ {Z 1 , Z 2 , Z 3 }, scenario in Table 2 and type of Copula, the larger τ , the smaller the control limit UCL Z . • For a fixed scenario in Table 2, value of τ and type of Copula, the upper control limits of the statistic Z ∈ {Z 1 , Z 2 , Z 3 } always satisfy UCL Z1 < UCL Z2 < UCL Z3 . • For a fixed scenario in Table 2, value of τ and statistic Z ∈ {Z 1 , Z 2 , Z 3 }, the values of UCL Z are more or less the same no matter the type of Copula considered.

e -θ -1 R\{0} τ = 1 + 4(D 1 (θ)-1)
As in Section 2.5 for independent times T and amplitudes X, EATS 1 values (for shifts in both T and X) have been reported in Tables 7-9 of [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Correlated Data[END_REF] in the case of dependent times T and amplitudes X. These Tables only consider the Frank's Copula and τ ∈ {0.2, 0.5, 0.8}. These results clearly show that, irrespective of the level of dependence, when a shift in both T and X is likeky to occur, the best option is to use the statistics Z 1 or, eventually, Z 2 , but not the statistic Z 3 which is never considered as an efficient one. These conclusions are similar to the ones obtained in Section 2.5 and they remain identical for other Copulas like the Clayton's or Gumbel's ones. More details concerning these aspects can be found in [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Correlated Data[END_REF].

Illustrative example

The following illustrative example has been detailed for the first time in [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Correlated Data[END_REF] and it is related to a company that has recorded for one of its bottleneck machine, during about 6 years (from 08/01/2012 to 27/12/2018) all the breakdown dates (D i in days) as well as the estimated corresponding incurred costs (X i , in euros) which include all the repair and restart costs (spare parts, manpower) and the cost of manufacturing disruption, see Table 8. This data set of 44 dates is divided into two parts

• 30 breakdowns recorded during 5 years (2012 to 2016) and used in this example as a Phase 1 data set. • 14 breakdowns recorded during 2 years (2017 to 2018) and used in this example as a Phase 2 data set.

The times T i and amplitude X i , i = 1, . . . , 44 in Table 8 have also been plotted in Figure 5 with • (for Phase 1) and • (for Phase 2), respectively. From the scatterplot shown in Figure 5, it can be noted that when the time T i between consecutive breakdowns is smaller (larger), the corresponding cost X i seems to be also smaller (larger), indicating a potential slight positive correlation between T and X. Investigations (during the period 2012 to 2016) about this phenomena have clarified why such a positive correlation between T and X may occur in this situation. After the occurrence of a breakdown, if the next one occurs after a short period of time, it is often due to the same problem occurring to the process: consequently, the time to looking for the breakdown causes is reduced. The spare parts costs are also reduced as they have already been purchased for the previous breakdown. On the contrary, when the next breakdown occurs after a long period of time, the causes are usually different from the previous breakdown and need more time to be revealed; new spare parts need to be purchased, thus increasing the cost. In order to evaluate if a positive correlation significantly exists between T and X, the Kendall's and Spearman's rank correlation coefficients have been estimated to τ = 0.4657 and ρ = 0.6129 as well as their correponding p-values 0.00035 and 0.00032, respectively, thus confirming a positive correlation between T and X. In this example, a Frank's Copula has been chosen to model the dependence between T and X. As explained in [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Correlated Data[END_REF], if τ = 0.4657 then an estimation of the dependence parameter is θ = 5.14.

Concerning the marginal distributions of T and X, the best fit using the Kolmogorov-Smirnov's test is to choose a Gamma distribution for the time T with parameters (a 0 = 11.6488, b 0 = 5.0562) and the Weibull distribution for the cost X with parameters (a 0 = 4.8472, b 0 = 5396.4958). 8. Assuming an in-control ATS value, ATS 0 = 9125 days (i.e. 25 years), the upper control limits of the three TBEA charts based on statistics Z 1 , Z 2 and Z 3 are found to be UCL Z1 = 0.57, UCL Z2 = 2.06 and UCL Z3 = 3.18, respectively, see [START_REF] Rahali | Evaluation of Shewhart Time-Between-Events-and-Amplitude Control Charts for Correlated Data[END_REF] for more details. The three TBEA charts, corresponding to the statistics Z 1 , Z 2 and Z 3 , are plotted in Figure 6 along with their upper control limits. As it can be seen, the Phase 1 part of these charts seems to confirm the fact that from 2012 to 2016, the time between consecutive breakdowns and their corresponding costs were in stable state. But, from 2017, things seems to have changed as several out-of-control situations (see values in bold in Table 8) have been detected by all the TBEA charts due to more frequent breakdowns and an increasing maintenance cost (also due to an aging machine). Every time an out-of-control situation is detected, the production has been stopped, the root causes of the breakdown have been searched for, analyzed and repaired. Then, the machine has been restarted.

A distribution-free TBEA chart

Motivation

The TBEA control charts developed in the previous sections are parametric ones, i.e. they assume that the distributions of the Times X and their Amplitudes X are perfectly known. However, as it is mentioned in [START_REF] Qiu | Introduction to Statistical Process Control[END_REF], in many practical situations, the distributions of these random variables are unknown (or their parameters cannot be reliably estimated by means of a Phase 1 retrospective study) making the implementation of traditional parametric control charts to be an incorrect approach. For instance, the parametric distributions that have been investigated in Sections 2 and 3, are the Gamma, Lognormal, Normal and Weibull, but are these distributions really appropriate in these cases? To overcome this problem, distribution-free control charts have been proposed and investigated in the literature. Among the most recent ones, we can cite for instance [START_REF] Celano | The Performance of the Shewhart Sign Control Chart for Finite Horizon Processes[END_REF], [START_REF] Abid | Use of Ranked Set Sampling in Nonparametric Control Charts[END_REF]Abid et al. ( , 2017aAbid et al. ( ,b, 2018)), [START_REF] Castagliola | An EWMA-type Sign Chart with Exact Run Length Properties[END_REF] and Chakraborti and Graham (2019b). Most of these approaches use nonparametric statistics like the Sign or the Wilcoxon signed-rank statistics. Practical guidelines for the implementation of such distribution-free control charts can be found in Chakraborti and Graham (2019a).

In this Section, we present an upper-sided distribution-free EWMA control chart for monitoring TBEA data. Moreover, as evaluating the Run Length properties of any EWMA scheme based on discrete data is a challenging problem, we will also introduce a specific method called "continuousify" which allows to obtain much better and replicable results.

Model

In this Section, we assume that F T (t|θ T ) and F X (x|θ X ) are the unknown continuous c.d.f. of T i and X i , i = 1, 2, . . ., where θ T and θ X are known quantiles, respectively. As for the previous sections, when the process is in-control, we have θ T = θ T0 , θ X = θ X0 and, when the process is out-of-control, we have θ T = θ T1 , θ X = θ X1 . Without loss of generality, we will consider that θ T and θ X are the median values of T i and X i , respectively. Of course, if necessary, 8.

other quantiles can also be considered.

Let

p T = P(T i > θ T0 |θ T ) = 1 -F T (θ T0 |θ T ) and p X = P(X i > θ X0 |θ X ) = 1 -F X (θ X0 |θ X ), i = 1, 2, .
. ., be the probabilities that T i and X i are larger than θ T0 and θ X0 assuming that the actual median values are θ T and θ X , respectively. Let us also define q T = 1p T and q X = 1p X .

If the process is in-control, we have p T = p T0 = 1 -F T (θ T0 |θ T0 ) = 0.5, p X = p X0 = 1 -F X (θ X0 |θ X0 ) = 0.5 and, when the process is out-of-control, we have p

T = p T1 = 1 -F T (θ T0 |θ T1 ), p X = p X1 = 1 -F X (θ X0 |θ X1 ).
The upper-sided distribution-free EWMA TBEA control chart that will be introduced in this section is based on the following sign statistics ST i and SX i , for i = 1, 2, . . . as

ST i = sign(T i -θ T0 ), SX i = sign(X i -θ X0 ),
where sign(x) = -1 if x < 0 and sign(x) = +1 if x > 0. The case x = 0 will not be considered here (even if it can happen in practice) because the random variables T i and X i are assumed to be continuous ones. In order to simultaneously monitor (T i , X i ), i = 1, 2, . . ., we suggest to define the statistic S i , for i = 1, 2, . . . as

S i = SX i -ST i 2 .
Since ST i ∈ {-1, +1} and SX i ∈ {-1, +1} we have S i ∈ {-1, 0, +1} and, more precisely, we have:

• S i = -1 when T i increases (ST i = +1
) and, at the same time, X i decreases (SX i = -1). In this case, the process seems to be in an acceptable situation. • S i = +1 when T i decreases (ST i = -1) and, at the same time, X i increases (SX i = +1). In this case, the process seems to be in an unacceptable situation.

• S i = 0 when both T i and X i increase or when both T i and X i decrease.

In this case, the process seems to be in an intermediate situation.

It can be easily proven that the p.m.f. (probability mass function) f Si (s|p T , p X ) = P(S i = s|p T , p X ) and the c.d.f. F Si (s|p T , p X ) = P(S i ≤ s|p T , p X ) of S i are equal to

f Si (s|p T , p X ) =        p T q X if s = -1 p T p X + q T q X if s = 0 q T p X if s = +1 0 if s ∈ {-1, 0, 1} , and
F Si (s|p T , p X ) =        0 if s ∈ (-∞, -1) p T q X if s ∈ [-1, 0) p T + q T q X if s ∈ [0, 1) 1 if s ∈ [1, +∞)
.

If we define an EWMA TBEA type control chart directly monitoring the statistic S i , it is known (see [START_REF] Wu | A Distribution-Free EWMA Control Chart for Monitoring Time-Between-Events-and-Amplitude Data[END_REF] for details) that, due to the strong discrete nature of this statistic, an accurate computation using Markov chain or integral equation methods, for instance, of the corresponding Run Length properties (ARL, SDRL, ...) is impossible. Consequently, it is actually possible to define an EWMA TBEA type control chart based on S i but it is unfortunately impossible to correctly evaluate its properties and, therefore, it is impossible to design it in a reliable way. In order to overcome this problem and before any implementation of an EWMA scheme, [START_REF] Wu | A Distribution-Free EWMA Control Chart for Monitoring Time-Between-Events-and-Amplitude Data[END_REF] suggested to define an extra parameter σ ∈ [0.1, 0.2] and to transform the discrete random variable S i into a new continuous one, denoted as S * i defined as a mixture of 3 normal random variables Y i,-1 ∼ Nor(-1, σ), Y i,0 ∼ Nor(0, σ) and Y i,+1 ∼ Nor(+1, σ), with weights w -1 = p T q X , w 0 = p T p X + q T q X and w +1 = q T p X (corresponding to the probabilities f Si (s|p T , p X ), s ∈ {-1, 0, +1}), respectively, i.e.

S * i =    Y i,-1 if S i = -1, Y i,0 if S i = 0, Y i,+1 if S i = +1.
This strategy has been called continuousify by [START_REF] Wu | A Distribution-Free EWMA Control Chart for Monitoring Time-Between-Events-and-Amplitude Data[END_REF]. It is easy to prove that the c.d.f.

F S * i (s|p T , p X ) = P(S * i ≤ s|p T , p X ) of S * i is equal to F S * i (s|p T , p X ) =p T q X F Nor (s| -1, σ) + (p T p X + q T q X )F Nor (s|0, σ) + q T p X F Nor (s| + 1, σ), (18)
and its expectation E(S * i ) and variance V(S * i ) are equal to

E(S * i ) = p X -p T , V(S * i ) = σ 2 + p T q T + p X q X .
When the process is in-control, we have p T0 = q T0 = 0.5, p X0 = q X0 = 0.5 and the expectation and variance of S * i simplify to E(S * i ) = 0 and V(S * i ) = σ 2 +0.5. As the main goal is to detect an increase in S i (in order to avoid, for instance, more damages or injuries/costs) rather than a decrease, the following uppersided EWMA TBEA control chart based on the statistic Z * i is proposed

Z * i = max(0, λS * i + (1 -λ)Z * i-1 ), (19)
with the following upper asymptotic control limit UCL defined as

UCL = E(S * i ) =0 +K λ 2 -λ × V(S * i ) = √ σ 2 +0.5 = K λ(σ 2 + 0.5) 2 -λ , (20)
where λ ∈ [0, 1] and K > 0 are the control chart parameters to be fixed and the initial value Z * 0 = 0. The zero-state ARL and SDRL of the proposed distribution-free upper-sided EWMA TBEA control chart can be obtained using the standard approach of [START_REF] Brook | An Approach to the Probability Distribution of CUSUM Run Length[END_REF] which assumes that the behavior of this control chart can be well represented by a discrete-time Markov chain with m + 2 states, where states i = 0, 1, . . . , m are transient and state m + 1 is an absorbing one. The transition probability matrix P of this discrete-time Markov chain is

P =   Q r 0 1   =        Q 0,0 Q 0,1 • • • Q 0,m r 0 Q 1,0 Q 1,1 • • • Q 1,m r 1 . . . . . . . . . . . . Q m,0 Q m,1 • • • Q m,m r m 0 0 • • • 0 1        ,
where Q is the (m + 1, m + 1) matrix of transient probabilities, where 0 = (0, 0, . . . , 0) and where the (m + 1, 1) vector r satisfies r = 1 -Q1 (i.e. row probabilities must sum to 1) with 1 = (1, 1, . . . , 1) . The transient states i = 1, . . . , m are obtained by dividing the interval [0, UCL] into m subintervals of width 2∆, where ∆ = UCL 2m . By definition, the midpoint of the i-th subinterval (representing state i) is equal to H i = (2i -1)∆. The transient state i = 0 corresponds to the "restart state" feature of our chart and it is represented by the value H 0 = 0. Concerning the proposed upper-sided EWMA TBEA control chart, it can be easily shown that the generic element Q i,j , i = 0, 1, . . . , m, of the matrix Q is equal to

• if j = 0, Q i,0 = F S * i - (1 -λ)H i λ p T , p X , (21)
• if j = 1, 2, . . . , m,

Qi,j = F S * i Hj + ∆ -(1 -λ)Hi λ pT , pX -F S * i Hj -∆ -(1 -λ)Hi λ pT , pX (22)
Let q = (q 0 , q 1 , . . . , q m ) be the (m + 1, 1) vector of initial probabilities associated with the m + 1 transient states. In our case, we assume q = (1, 0, . . . , 0) , i.e. the initial state corresponds to the "restart state". When the number m of subintervals is sufficiently large (say m = 300), this finite approach provides an effective method that allows the ARL and SDRL to be accurately evaluated using the following classical formulas

ARL = q (I -Q) -1 1, ( 23
) SDRL = 2q (I -Q) -2 Q1 + ARL(1 -ARL). ( 24
)
In order to illustrate the "power" of the continuousify technique on the uppersided EWMA TBEA control chart (in the case K = 3 and λ = 0.2), Table 9 presents some ARL values obtained without (left side) and with (right side) this technique, corresponding to 3 combinations for (p X , p T ), a number of subintervals m ∈ {100, 120, . . . , 400} and σ = 0.125. When the continuousify technique is not used, the random variable S * i in ( 19) is replaced by S i or, equivalently, the parameter σ is set to 0. For comparison purpose, the last row of Table 9 also provides the ARL values obtained by simulations. As it can be seen, if the continuousify technique is not used (left side), the ARL values obtained using the Markov chain method can have a large variability without any clear monotonic convergence when m increases. It turns out that with these unstable ARL values, it can be really difficult to design and optimize the upper-sided EWMA TBEA control chart. On the other side, if the continuousify technique is used (right side), the ARL values obtained using the Markov chain method are clearly very stable, even for small values of m, but (and this is the price to pay for this stability) they are a bit larger than those obtained by simulations (compare 26.08, 12.23, 27.88 vs. 24.71, 11.66, 26.46). This property is due to the extra term σ > 0 in (20).

Comparative studies

Since the continuousify technique is demonstrated to provide reliable ARL values, it is therefore possible to compute the optimal chart parameters (λ * , K * ) for the upper-sided EWMA TBEA control chart minimizing the out-of-control ARL(λ * , K * , σ, p T , p X ) for p T = 0.5 and p X = 0.5 under the constraint ARL(λ * , K * , σ, 0.5, 0.5) = ARL 0 , where ARL 0 is a predefined value for the in-control ARL. These optimal values are listed in Table 10 with the corresponding out-of-control (ARL, SDRL) values for p T ∈ {0.1, 0.2, . . . , 0.4} (only considering a decrease in T ), p X ∈ {0.5, 0.6, . . . , 0.9} (only considering an increase in X), for four possible choices for σ ∈ {0.1, 0.125, 0.15, 0.2} and assuming ARL 0 = 370.4. These values of (λ * , K * ) can be freely be used by practitioners who need to optimally detect a specific shift in the times and/or in the amplitudes.

A comparison between the upper-sided EWMA TBEA control chart introduced in this Section and the three parametric TBEA control charts presented in Subsection 2.1 has been investigated in [START_REF] Wu | A Distribution-Free EWMA Control Chart for Monitoring Time-Between-Events-and-Amplitude Data[END_REF] using the EARL 1 (instead of the EATS 1 ) for the following two scenarios

• Scenario #1: a Normal distribution for X with in-control mean µ X0 = 10 and standard-deviation σ X0 = 1 and a gamma distribution for T with incontrol mean µ T0 = 10 and standard-deviation σ T0 = 2, i.e. X ∼ Nor(10, 1) and T ∼ Gam(25, 0.4). • Scenario #2: a Normal distribution for X with in-control mean µ X0 = 10 and standard-deviation σ X0 = 2 and a Weibull distribution for T with incontrol mean µ T0 = 10 and standard-deviation σ T0 = 1, i.e. X ∼ Nor(10, 2) and T ∼ Wei (12.1534, 10.4304).

Based on the results in Tables 3 and4 in [START_REF] Wu | A Distribution-Free EWMA Control Chart for Monitoring Time-Between-Events-and-Amplitude Data[END_REF], the conclusion is that no matter the scenario (#1 or #2) or the statistic considered Z ∈ {Z 1 , Z 2 , Z 3 }, the out-of-control EARL 1 values obtained for the distributionfree upper-sided EWMA TBEA control chart are always smaller than the ones obtained for the parametric Shewhart control charts introduced in Subsection 2.1, thus showing the advantage of using the proposed distribution-free control chart in situations where the distributions for T and X are unknown.

Illustrative example

We consider here the same illustrative example as the one already presented in Section 2.6 concerning the days T i between fires in forests of the french region "Provence -Alpes -Côte D'Azur" and their amplitudes X i (burned surface in ha = 10000m 2 ). In order to compute the control limit UCL of the distribution-free upper-sided EWMA TBEA chart, the following values have been fixed: p T = 0.3, p X = 0.7, σ = 0.125 and ARL 0 = 370.4. The corresponding optimal values for λ and K are found to be λ * = 0.07 and The in-control median values for T and X have been estimated from the Phase 1 data set and they are equal to θ T0 = 3 days and θ X0 = 5.3 ha. These values are used to compute the values ST i , SX i , S i and S * i in Table 11. As it can be noticed, some dates are such that T iθ T0 = 0. Of course, this situation is not supposed to happen as the times T i are supposed to be continuous random variables but, due to the measurement scale (days), this situation actually happens. When this situation occurs, a possible simple strategy consists in assigning ST i = 0 (instead of -1 or +1). For this reason, in Table 11, some values of S i = s = ±0.5 and the corresponding values for S * i are obtained by randomly generating a Nor(s, σ) random variable, as it is already the case for values s ∈ {-1, 0, +1}. For instance, in Table 11, when D i = 70 we have S i = 0.5 and the corresponding value for S * i has been randomly generated from a Nor(0.5, 0.125) distribution (S * i = 0.552). The values Z * i have been computed using eqn.( 19), for both Phase 1 and 2 data sets, recorded in Table 11 and plotted in Figure 7 along with the distribution-free EWMA TBEA upper control limit UCL = 0.344. If the distribution-free upper-sided EWMA TBEA chart does not detect any out-of-control situations during the Phase 1 (validating the in-control state of this phase), it nevertheless detects several out-of-control situations during the period mid-June 2017 -end of September 2017, (see also the bold values in Table 11), confirming that a decrease in the time between fires occurred with a concurrent increase in the amplitude of these fires. This conclusion is consistent with the one obtained in Section 2.6 in which a parametric approach assuming a lognormal distribution for both T i and X i was used.

Conclusions

The three contributive Sections of this Chapter have clearly demonstrated that efficient solutions do exist when the aim is to simultaneously monitor the time interval T of an event E as well as its amplitude X. These solutions can be either parametric, for independent or dependent situations, and they can also be distribution-free if there is no a priori knowledge about the distributions associated with T and X.

In our opinion, future researches on the monitoring of TBEA data can be undertaken toward several directions:

• In the proposed parametric approaches, the estimation of the parameters (for the distributions or the Copulas) is not taken into account at all in the design and evaluation of the TBEA control charts. The impact of the Di Ti Xi STi SXi Si S * i Z * i 9 9 3.68 1 -1 -1.0 -0.917 0.000 258 1 1.00 -1 -1 0.0 -0.078 0.000 26 17 1.99 1 -1 -1.0 -0.802 0.000 260 2 3.70 -1 -1 0.0 0.119 0.008 60 34 6.00 1 1 0.0 -0.081 0.000 262 2 3.17 -1 -1 0.0 -0.063 0.003 67 7 1.19 1 -1 -1.0 -0.901 0.000 265 3 18.40 0 1 0.5 0.333 0.026 70 3 135.80 0 1 0.5 0.552 0.039 268 3 1.00 0 -1 -0.5 -0.145 0.014 72 2 14.37 -1 1 1.0 1.113 0.114 269 1 2.22 -1 -1 0.0 0.208 0.028 86 14 8.10 1 1 0.0 -0.104 0.099 271 2 19.09 -1 1 1.0 1.001 0.096 88 2 32.31 -1 1 1.0 0.892 0.154 272 1 2.00 -1 -1 0.0 0.027 0.091 94 6 3.07 1 -1 -1.0 -1.056 0.069 274 2 34.28 -1 1 1.0 1.086 0.161 95 1 10.03 -1 1 1.0 0.867 0.125 276 2 3.00 -1 -1 0.0 0.070 0.154 96 1 7.93 -1 1 1.0 1.033 0.189 277 1 6.63 -1 1 1.0 0.955 0.210 97 1 1.50 -1 -1 0.0 0.409 0.204 278 1 4.47 -1 -1 0.0 -0.097 0.189 103 6 23.30 1 1 0.0 -0.116 0.182 285 7 8.24 1 1 0.0 0.160 0.187 106 3 3.73 0 -1 -0.5 -0.708 0.120 286 1 769.45 -1 1 1.0 1.024 0.246 109 3 4.73 0 -1 -0.5 -0.677 0.064 287 1 4.37 -1 -1 0.0 -0.144 0.218 111 2 3.19 -1 -1 0.0 0.179 0.072 288 1 90.70 -1 1 1.0 0.961 0.270 113 2 6.25 -1 1 1.0 1.032 0.139 289 1 11.49 -1 1 1.0 1.044 0.324 114 1 3.60 -1 -1 0.0 -0.155 0.118 295 6 3590.78 1 1 0.0 0.033 0.304 115 1 6.12 -1 1 1.0 1.112 0.188 296 1 1427.92 -1 1 1.0 0.949 0.349 118 3 1.50 0 -1 -0.5 -0.740 0.123 297 1 255.96 -1 1 1.0 1.054 0.399 122 4 1.33 1 -1 -1.0 -1.009 0.044 298 1 1.00 -1 -1 0.0 -0.051 0.367 134 12 1.42 1 -1 -1.0 -1.037 0.000 302 4 13.88 1 1 0.0 -0.074 0.336 137 3 5.75 0 1 0.5 0.629 0.044 303 1 138.28 -1 1 1.0 1.117 0.391 140 3 3.47 0 -1 -0.5 -0.507 0.005 305 2 8.90 -1 1 1.0 1.153 0.444 142 2 13.31 -1 1 1.0 1.217 0.090 308 3 1.50 0 -1 -0.5 -0.342 0.389 143 1 26.31 -1 1 1.0 1.041 0.157 312 4 34.63 1 1 0.0 -0.217 0.347 144 1 18.54 -1 1 1.0 0.923 0.210 313 1 82.56 -1 1 1.0 0.811 0.379 146 2 66.17 -1 1 1.0 1.124 0.274 314 1 2.00 -1 -1 0.0 -0.019 0.351 147 1 9.90 -1 1 1.0 0.916 0.319 315 1 162.08 -1 1 1.0 1.071 0.402 150 3 4.22 0 -1 -0.5 -0.534 0.260 319 4 3.26 1 -1 -1.0 -1.056 0.300 157 7 34.28 1 1 0.0 -0.110 0.234 321 2 285.91 -1 1 1.0 0.729 0.330 161 4 2.23 1 -1 -1.0 -1.102 0.140 322 1 2.00 -1 -1 0.0 -0.283 0.287 162 1 1.84 -1 -1 0.0 0.152 0.141 325 3 11.57 0 1 0.5 0.347 0.291 163 1 2.88 -1 -1 0.0 -0.018 0.130 334 9 34.70 1 1 0.0 0.068 0.275 164 1 21.46 -1 1 1.0 1.087 0.197 335 1 431.00 -1 1 1.0 1.150 0.337 165 1 4.46 -1 -1 0.0 -0.001 0.183 336 1 10.89 -1 1 1.0 1.003 0.383 166 1 58.27 -1 1 1.0 1.034 0.243 340 4 1.00 1 -1 -1.0 -1.004 0.286 167 1 8.84 -1 1 1.0 0.863 0.286 346 6 1.50 1 -1 -1.0 -0.921 0.202 180 13 1.03 1 -1 -1.0 -0.905 0.203 347 1 1.17 -1 -1 0.0 0.100 0.195 187 7 16.57 1 1 0.0 0.156 0.199 349 2 1.27 -1 -1 0.0 0.129 0.190 201 14 4.96 1 -1 -1.0 -1.084 0.110 350 1 26.25 -1 1 1.0 1.098 0.254 202 1 1.37 -1 -1 0.0 -0.087 0.096 353 3 11.66 0 1 0.5 0.332 0.259 205 3 23.39 0 1 0.5 0.498 0.124 354 1 3.03 -1 -1 0.0 0.127 0.250 225 20 1.70 1 -1 -1.0 -1.032 0.043 355 1 12.00 -1 1 1.0 1.130 0.311 247 22 5.30 1 0 -0.5 -0.727 0.000 356 1 1.00 -1 -1 0.0 -0.206 0.275 248 1 15.64 -1 1 1.0 1.161 0.081 257 9 5.14 1 -1 -1.0 -0.921 0.011 11 parameter estimation is known to strongly influence the efficiency of any control chart and, therefore, researches on this topic should be done.

• Measures like times or amplitudes are obviously subject to measurement errors. These kinds of error are also known to negatively impact the efficiency of any control chart. Reasearches investigating the impact of the measurement errors on T and/or X should also be undertaken in order to evaluate how much they actually impact the performance of parametric TBEA control charts. • Instead of considering X as univariate random variable, it could be considered in some cases as a p-variate random vector X = (X 1 , . . . , X p ) where each X k is the amplitude of a specific characteristic and the goal would be to simultaneously monitor (T, X). For instance, in the forest fires example, the amplitude could be considered as a bivariate vector X = (X 1 , X 2 ) where X 1 would be the burned surface and X 2 would be the cost related to the fires. • Often, historical data availability in monitoring of adverse events is limited to a few records. Thus, knowledge about the frequency distribution of these events is too restricted to fit a reliable statistical model. With these scenarios, there is room for approaching the monitoring problem with distribution-free approaches, which, therefore, deserve a lot of attention by researchers.

Fig. 1 .

 1 Fig. 1. Dates of occurrence D0 = 0, D1, D2, . . ., time intervals T1 = D1 -D0, T2 = D2 -D1, . . . and amplitudes X1, X2, . . . of a negative event E

Fig. 2 .

 2 Fig. 2. Integration areas used for statistics (a) and (b) Z1, (c) Z2 and (d) Z3

θ

  Note: D1(θ) is the Debye function of the first kind defined as D1(θ)

Fig. 5 .

 5 Fig.5. Phase 1 (•) and 2 (•) data corresponding to the time (T in days) between breakdowns and amplitudes (X in euros) corresponding to the data set in Table8.

Fig. 6 .

 6 Fig.6. Statistics Z1, Z2 and Z3 corresponding to the data set in Table8.

  without continuousify with continuousify (σ = 0.125) pT = 0.3 pT = 0.2 pT = 0.1 pT = 0.3 pT = 0.2 pT = 0.1 m pX = 0.8 pX = 0.9 pX = 0.6 pX = 0.8 pX = 0.9 pX = 0.

Fig. 7 .

 7 Fig. 7. Distribution-free EWMA TBEA chart with statistic Z * i corresponding to the data set in Table 11

C

  .W. Zhang, M. Xie, J.Y. Liu, and T.N. Goh. A Control Chart for the Gamma Distribution as a Model of Time Between Events. International Journal of Production Research, 45(23):5649-5666, 2007. H.Y. Zhang, M. Shamsuzzaman, M. Xie, and T.N. Goh. Design and Application of Exponential Chart for Monitoring Time-Between-Events Data under Random Process Shift. The International Journal of Advanced Manufacturing Technology, 57(9):849-857, 2011a. H.Y. Zhang, M. Xie, T.N. Goh, and M. Shamsuzzaman. Economic Design of Time-Between-Events Control Chart System. Computers & Industrial Engineering, 60(4):485-492, 2011b.

Table 1 .

 1 Distributions used for the comparison of the 3 TBEA charts

	Names	Parameters	f (x|a, b)

Table 2 .

 2 In-control configurations to be investigated

	Distributions T X	a0	b0	µ0 σ0 γ0
	Gamma	• • • • • •	100 25 4	0.1 10 1 0.4 10 2 2.5 10 5	0.2 0.4 1
	Lognormal • • -23.0334 10.0249 10 1 0.3010 • • -11.5277 5.0494 10 2 0.6080 • • -4.6382 2.1169 10 5 1.6250
	Normal	• • • •	10 10	1 2	10 1 10 2	0 0
	Weibull	• • 12.1534 10.4304 10 1 -0.7155 • • 5.7974 10.7998 10 2 -0.3519 • • 2.1013 11.2906 10 5 0.5664

Table 3 .

 3 Upper control limits UCLZ for the three TBEA charts based on statistics Z ∈ {Z1, Z2, Z3}, satisfying ATS0 = 370.4

	Statistic Z1

Table 5 .

 5 Phase 1 and 2 data sets corresponding to time (T in days) between fires, amplitudes (X as the burned surface in ha) and the values of the statistics Z1, Z2 and Z3.

			Phase 1				Phase 2
	i T	X	Z1 Z2	Z3	i T	X	Z1	Z2	Z3
	1 3.68 -1.37 0.16 0.88	1 1 1.00	-0.11	0.40	5.54
	2 17 1.99 -2.96 0.05 0.47	2 2 3.70	-0.09	0.75	3.01
	3 34 6.00 -5.78 0.07 0.60	3 2 3.17	-0.13	0.64	2.97
	4 1.19 -1.19 0.07 0.87	4 3 18.40	0.81	2.47	3.18
	5 135.80 9.45 18.23 11.82	5 3 1.00	-0.47	0.13	1.90
	6 14.37 0.69 2.89 3.79	6 1 2.22	-0.02	0.89	5.63
	7 14 8.10 -1.96 0.23 0.99	7 2 19.09	1.04	3.84	4.14
	8 32.31 2.01 6.51 5.11	8 1 2.00	-0.04	0.81	5.62
	9 3.07 -0.87 0.21 1.14	9 2 34.28	2.16	6.90	5.26
	10 10.03 0.56 4.04 6.21	10 2 3.00	-0.14	0.60	2.95
	11 7.93 0.40 3.19 6.05	11 1 6.63	0.31	2.67	5.96
	12 1.50 -0.07 0.60 5.58	12 1 4.47	0.15	1.80	5.80
	13 23.30 0.62 1.56 2.63	13 7 8.24	-0.67	0.47	1.39
	14 3.73 -0.27 0.50 2.10	14 1 769.45 56.49 309.87 62.14
	15 4.73 -0.20 0.63 2.17	15 1 4.37	0.14	1.76	5.79
	16 3.19 -0.13 0.64 2.97	16 1 90.70 6.50 36.53 12.15
	17 6.25 0.09 1.26 3.19	17 1 11.49	0.66	4.63	6.31
	18 3.60 0.08 1.45 5.73	18 6 3590.78 263.36 241.01 265.37
	19 6.12 0.27 2.46 5.92	19 1 1427.92 104.98 575.04 110.63
	20 1.50 -0.44 0.20 1.93	20 1 255.96 18.67 103.08 24.32
	21 1.33 -0.63 0.13 1.46	21 1 1.00	-0.11	0.40	5.54
	22 12 1.42 -2.09 0.05 0.56	22 4 13.88	0.29	1.40	2.39
	23 5.75 -0.13 0.77 2.25	23 1 138.28 10.00 55.69 15.65
	24 3.47 -0.29 0.47 2.08	24 2 8.90	0.29	1.79	3.39
	25 13.31 0.61 2.68 3.71	25 3 1.50	-0.44	0.20	1.93
	26 26.31 1.75 10.60 7.41	26 4 34.63	1.82	3.49	3.92
	27 18.54 1.18 7.47 6.83	27 1 82.56	5.90 33.25 11.55
	28 66.17 4.51 13.32 7.61	28 1 2.00	-0.04	0.81	5.62
	29 9.90 0.55 3.99 6.20	29 1 162.08 11.75 65.27 17.40
	30 4.22 -0.24 0.57 2.13	30 4 3.26	-0.49	0.33	1.61
	31 34.28 1.24 1.97 3.31	31 2 285.91 20.69 57.57 23.79
	32 2.23 -0.57 0.22 1.53	32 1 2.00	-0.04	0.81	5.62
	33 1.84 -0.05 0.74 5.60	33 3 11.57	0.30	1.55	2.67
	34 2.88 0.03 1.16 5.68	34 9 34.70	0.91	1.55	3.16
	35 21.46 1.40 8.64 7.05	35 1 431.00 31.56 173.57 37.21
	36 4.46 0.15 1.80 5.80	36 1 10.89	0.62	4.39	6.27
	37 58.27 4.11 23.47 9.76	37 4 1.00	-0.66	0.10	1.44
	38 8.84 0.47 3.56 6.12	38 6 1.50	-0.99	0.10	1.02
	39 13 1.03 -2.30 0.03 0.50	39 1 1.17	-0.10	0.47	5.55
	40 16.57 -0.06 0.95 2.00	40 2 1.27	-0.27	0.26	2.83
	41 14 4.96 -2.20 0.14 0.76	41 1 26.25	1.75	10.57	7.40
	42 1.37 -0.08 0.55 5.57	42 3 11.66	0.31	1.57	2.68
	43 23.39 1.17 3.14 3.55	43 1 3.03	0.04	1.22	5.69
	44 20 1.70 -3.53 0.03 0.40	44 1 12.00	0.70	4.83	6.35
	45 22 5.30 -3.63 0.10 0.64	45 1 1.00	-0.11	0.40	5.54
	46 15.64 0.97 6.30 6.62
	47 5.14 -1.27 0.23 0.99

  50

					Phase 1							Phase 2
	1
	Z
												U CL
	Oct	Nov	Dec	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep
	16	16	16	17	17	17	17	17	17	17	17	17
					Phase 1							Phase 2
	2
	Z
												U CL
	Oct	Nov	Dec	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep
	16	16	16	17	17	17	17	17	17	17	17	17
					Phase 1							Phase 2
	Z
												U CL
	Oct	Nov	Dec	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep
	16	16	16	17	17	17	17	17	17	17	17	17

3 Fig. 4. Statistics Z1, Z2 and Z3 corresponding to the data set in Table 5

Table 6 .

 6 Details concerning Gumbel's,Clayton's and Frank's Copulas

	Name	C(u, v|θ)	Domain for θ	τ and θ
	Gumbel exp -(-ln(u)

Table 7 .

 7 Pre-computed values of θ for several selected values of τ ∈

	{0, 0.1, 0.2, . . . , 0.9}
		θ
	τ Frank Clayton Gumbel
	0.0 0.00	0.00	1.00
	0.1 0.91	0.22	1.11
	0.2 1.86	0.50	1.25
	0.3 2.92	0.86	1.43
	0.4 4.16	1.33	1.67
	0.5 5.74	2.00	2.00
	0.6 7.93	3.00	2.50
	0.7 11.41 4.67	3.33
	0.8 18.19 8.00	5.00
	0.9 38.28 18.00	10.00

Table 8 .

 8 Phase 1 and 2 data sets corresponding to the time (Ti in days) between two consecutive breakdowns, amplitudes (Xi cost in euros) and the values of the statistics Z1, Z2 and Z3.

		Phase 1		Phase 2
	Date	i Ti Xi Z1,i Z2,i Z3,i	Date	i Ti Xi Z1,i Z2,i Z3,i
	10/03/12 1 4890 -0.064 0.939 1.939	11/01/17 1 63 5080 -0.043 0.960 1.962
	28/05/12 2 6180 -0.092 0.932 1.995	21/03/17 2 69 5350 -0.090 0.923 1.935
	25/07/12 3 3730 -0.231 0.766 1.770	07/05/17 3 47 3770 -0.036 0.955 2.015
	27/08/12 4 2930 0.032 1.057 2.377	15/07/17 4 69 4590 -0.243 0.792 1.782
	20/11/12 5 7600 0.093 1.065 2.230	14/10/17 5 91 5940 -0.344 0.777 1.848
	20/02/13 6 5580 -0.434 0.722 1.768	18/12/17 6 65 5420 -0.008 0.993 2.002
	30/04/13 7 4570 -0.247 0.789 1.778	26/02/18 7 70 4580 -0.262 0.779 1.767
	06/07/13 8 5230 -0.080 0.930 1.937	21/04/18 8 54 5430 0.181 1.197 2.189
	18/08/13 9 4470 0.174 1.238 2.274	14/05/18 9 23 5740 0.770 2.972 3.721
	28/09/13 10 3420 -0.005 0.993 2.128	27/06/18 10 44 6110 0.488 1.654 2.574
	22/11/13 11 3460 -0.234 0.749 1.770	22/08/18 11 56 7340 0.533 1.561 2.536
	08/02/14 12 5360 -0.241 0.818 1.839	29/10/18 12 68 8160 0.495 1.429 2.516
	05/04/14 13 4470 -0.047 0.951 1.956	24/11/18 13 26 4800 0.529 2.199 3.236
	28/05/14 14 4470 0.004 1.004 2.015	27/12/18 14 33 6570 0.768 2.371 3.113
	08/07/14 15 3320 -0.025 0.964 2.108
	27/09/14 16 4910 -0.382 0.722 1.720
	29/10/14 17 5010 0.470 1.864 2.854
	07/01/15 18 6630 0.152 1.128 2.182
	30/03/15 19 5710 -0.238 0.829 1.873
	20/05/15 20 5130 0.171 1.198 2.192
	16/07/15 21 5330 0.110 1.114 2.111
	01/09/15 22 5010 0.215 1.269 2.266
	22/10/15 23 3660 -0.126 0.855 1.895
	15/11/15 24 3340 0.268 1.657 3.129
	12/01/16 25 3600 -0.257 0.739 1.743
	14/03/16 26 5560 0.072 1.068 2.074
	28/04/16 27 5760 0.401 1.524 2.473
	25/06/16 28 6440 0.317 1.322 2.318
	16/08/16 29 6310 0.393 1.445 2.408
	09/11/16 30 6300 -0.169 0.883 1.967

Table 9 .

 9 ARL for the distribution-free EWMA TBEA chart computed with and without the continuousify technique.

Table 10 .

 10 Optimal values for (λ * , K * ) with the corresponding out-of-control values of (ARL, SDRL) for pT ∈ {0.1, 0.2, . . . , 0.4}, pX ∈ {0.5, 0.6, . . . , 0.9} and σ ∈ {0.1, 0.125, 0.15, 0.2}

				σ = 0.1
				pX
	pT	0.5	0.6	0.7	0.8	0.9
	0.5	(-,-)
		(370.40,-)
	0.4 (0.010,1.773) (0.025,2.174)
		(105.66,74.04) (50.77,32.32)
	0.3 (0.025,2.174) (0.045,2.387) (0.070,2.515)
		(51.54,32.55) (30.55,18.04) (20.50,11.38)
	0.2 (0.040,2.348) (0.070,2.515) (0.100,2.591) (0.145,2.639)
		(31.30,17.51) (20.74,11.40) (14.85,7.67) (11.19,5.55)
	0.1 (0.060,2.474) (0.090,2.571) (0.135,2.634) (0.180,2.645) (0.240,2.627)
		(21.40,10.76) (15.16,7.37) (11.32,5.40) (8.76,3.84) (6.99,2.74)
				σ = 0.125
				pX
	pT	0.5	0.6	0.7	0.8	0.9
	0.5	(-,-)
		(370.40,-)
	0.4 (0.010,1.774) (0.025,2.174)
		(106.19,74.55) (51.11,32.63)
	0.3 (0.025,2.174) (0.045,2.387) (0.070,2.515)
		(51.88,32.87) (30.79,18.25) (20.68,11.53)
	0.2 (0.040,2.348) (0.065,2.496) (0.100,2.592) (0.140,2.638)
		(31.53,17.72) (20.91,11.27) (14.99,7.80) (11.32,5.57)
	0.1 (0.060,2.474) (0.090,2.572) (0.135,2.634) (0.175,2.648) (0.225,2.639)
		(21.57,10.92) (15.30,7.50) (11.44,5.49) (8.88,3.89) (7.10,2.75)
				σ = 0.15
				pX
	pT	0.5	0.6	0.7	0.8	0.9
	0.5	(-,-)
		(370.40,-)
	0.4 (0.010,1.775) (0.025,2.175)
		(106.83,75.16) (51.53,33.01)
	0.3 (0.025,2.175) (0.045,2.387) (0.070,2.515)
		(52.30,33.26) (31.08,18.51) (20.90,11.71)
	0.2 (0.040,2.348) (0.065,2.496) (0.095,2.584) (0.135,2.636)
		(31.82,17.97) (21.13,11.45) (15.17,7.81) (11.47,5.61)
	0.1 (0.055,2.449) (0.090,2.573) (0.130,2.632) (0.170,2.651) (0.215,2.646)
		(21.79,10.76) (15.47,7.64) (11.59,5.53) (9.02,3.96) (7.23,2.80)
				σ = 0.2
				pX
	pT	0.5	0.6	0.7	0.8	0.9
	0.5	(-,-)
		(370.40,-)
	0.4 (0.010,1.777) (0.025,2.176)
		(108.43,76.68) (52.57,33.96)
	0.3 (0.020,2.085) (0.045,2.387) (0.065,2.496)
		(53.33,32.51) (31.81,19.15) (21.44,11.90)
	0.2 (0.040,2.348) (0.065,2.496) (0.090,2.574) (0.125,2.630)
		(32.53,18.61) (21.66,11.92) (15.60,8.01) (11.84,5.74)
	0.1 (0.055,2.449) (0.085,2.562) (0.120,2.624) (0.155,2.652) (0.195,2.658)
		(22.31,11.21) (15.89,7.83) (11.95,5.64) (9.34,4.06) (7.53,2.91)

Table 11 .

 11 Phase 1 and 2 values of Di, Ti, Xi, STi, SXi, Si, S * i and Z * i for the forest fires example.

	Phase 1	Phase 2
	Di Ti Xi STi SXi Si S * i	Z * i

Philippe Castagliola, Giovanni Celano, Dorra Rahali, and Shu Wu