César, Mohand Je Voudrais Aussi Remercier Ansaar

La méthode atteint également un F1-score supérieur à 99% avec des mesures terrain représentant cinq classes sur trois sites industriels différents. Enfin, nous étudions également l'impact des creux de tension sur les équipements industriels. Nous proposons une méthodologie pour estimer la composition de charges déconnectées suite à un creux de tension. Les résultats ont montré des limites en termes d'interaction harmonique entre les charges. Ces limites sont discutées et plusieurs propositions sont faites pour améliorer la méthode.

À l'issue de ce travail de thèse, je souhaiterais remercier tous ceux qui auront contribué à la réussite de ces trois années de recherche.

Je tiens tout d'abord à remercier mon directeur de thèse Claude Delpha, qui m'a encadré et suivi tout au long de cette thèse, pour sa grande gentillesse, sa disponibilité permanente et pour ses conseils précieux aux moments charniers. Je tiens également à remercier très chaleureusement mes encadrants Demba Diallo, Sophie Bercu et Ludovic Bertin, pour m'avoir guidé tout au long de mes travaux. Merci de m'avoir accordé votre temps, de m'avoir enseigné la rigueur, tout en restant dans la bienveillance et toujours avec le sourire. J'adresse également mes remerciements à M. Guy Clerc et M. Ghaleb Hoblos, pour l'honneur qu'ils m'ont fait en acceptant d'être rapporteurs de cette thèse et pour leurs remarques très pertinentes, ainsi qu'à tous les membres du jury pour la qualité de leurs questions et de nos échanges lors de ma soutenance.

Ce travail n'aurait pas été possible sans le soutien financier et industriel d'EDF R&D. Je tiens ainsi à remercier particulièrement mon manager Mathieu Caujolle, qui m'a donné l'opportunité de réaliser cette thèse au sein de son équipe, et m'a accompagné depuis les premiers jours de ce projet. Merci de m'avoir accordé ta confiance et pour tes nombreux encouragements, même pendant les moments les plus difficiles.

Titre : Automatisation du diagnostic de la qualité d'électricité des réseaux électriques industriels Mots clés : Qualité d'électricité, réseaux électriques industriels, creux de tension, classification de séries temporelles multivariées Résumé : L'analyse de la qualité de l'électricité est une demande qui s'est vue augmentée au cours des dernières décennies. Les creux de tension sont les perturbations les plus fréquentes et les plus impactantes dans les réseaux électriques industriels, entraînant des pertes financières importantes pour les clients industriels. Le coeur de ce travail de thèse est dédié à la classification des causes de creux de tension ainsi qu'à leur localisation relative par rapport au point de mesure principale. L'algorithme développé utilise les formes d'onde de tension et de courant comme entrées pour identifier les causes des creux de tension dans les réseaux industriels BT. La solution est basée sur des signatures de séries temporelles quadridimensionnelles, obtenues par l'application de la transformée de Fourier court terme (STFT) et de la transformée de Fortescue. La source d'un creux de tension est identifiée à l'aide d'une stratégie de classification basée sur la distance avec une mesure basée sur l'algorithme Dynamic Time Warping (DTW). En outre, l'algorithme soft-DTW est utilisé pour réduire la taille de la base de signatures d'apprentissage, augmentant ainsi la vitesse de classification. Les performances de la méthode ont été analysées en termes de séparabilité des classes, d'efficacité de la prédiction (précision et robustesse au bruit) et de sensibilité aux variations de la fréquence fon-Title : Automation of power quality diagnosis of industrial electrical grids Keywords : Power quality, industrial grids, voltage sag, multivariate time series classification Abstract : The demand for power quality analysis has increased over the past decades. Voltage sags are the most frequent and impactful disturbances in industrial power grids, leading to high financial losses for industrial clients. The core of this thesis work is dedicated to the classification of voltage sag causes and their relative location to the monitoring point. The solution uses voltage and current waveforms as input to identify the causes of voltage sags in LV industrial grids. The methodology is based on four-dimension time series signatures, obtained through the application of the Short-Time Fourier Transform (STFT) and the Fortescue Transform. The source of a voltage sag is identified using a distance-based classification strategy with a custom distance measure based on the Dynamic Time Warping algorithm (DTW). In addition, the soft-DTW algorithm is used to reduce the size of the signature training database and increase speed. The performance of the method was analyzed in terms of class separability, prediction efficiency (accuracy and robustness to noise), and sensitivity to fundamental frequency variations. The proposal is resilient regarding noise levels up to an SNR = 15 dB and fundamental frequency variations up to a shifting value ϵ = ±0.5Hz. Moreover, a confidence index on the prediction is proposed, increasing the algorithm's reliability. The proposal offers an easy implementation in industrial applications with no previous recorded data. It has the benefit of using a reduced-size reference database, entirely composed of synthetic data. The main advantages of the proposed method are its generalization capabilities and the possibility of raising an alert based on the confidence index. The obtained classification accuracy on synthetic data with seven causes is 100%. The method reaches a classification F1-score higher than 99% with field measurements representing five classes obtained from three different industrial sites. Finally, we also study the impact of voltage sags on industrial equipment. We propose a methodology to estimate the self-disconnected load composition following a voltage sag. The results showed some limitations in terms of harmonic interaction among the loads. Some of the limits of this approach are discussed and several proposals to improve the load composition estimation for future work are made.

General Conclusion Résumé en français

List of publications

International conference papers

Background and Motivations

The increase in electricity consumption, new usages, and the need to preserve the environment are at the origin of the energy transition. This transition is profoundly modifying electrical networks with the greater penetration of renewable energies associated with power electronics interfaces in the energy mix and control digitalisation [START_REF] Steimer | Enabled by high power electronics -Energy efficiency, renewables and smart grids[END_REF][START_REF] Tang | Power electronics: The enabling technology for renewable energy integration[END_REF]. Energy flows are becoming multidirectional and require more intelligent control to meet the electrical power supply's availability, reliability, safety and quality.

Power quality analysis has become an increasing concern to energy suppliers and their customers in the last few years. A reliable energy supply ensures regular operation of the electrical equipment in the network. Poor power quality induce the disruption of production lines or services, equipment malfunction, or even equipment damage, resulting in important financial losses [START_REF] Targosz | Pan-European power quality survey[END_REF] for industrial and large tertiary customers such as hospitals, data centers, etc. Indeed, energy quality is an important issue for manufacturers facing global competition. Thus, identifying the origin of power quality disturbances and assessing their effects on industrial equipment is essential for finding and providing adapted and cost-effective mitigation solutions to reduce the impact on the site's productivity.

Today, the analysis of power quality disturbances is most usually performed by experts in the field. It requires a high level of knowledge and expertise to make a reliable diagnosis, and propose relevant solutions. However, this process is highly time-consuming and requires specific intervention on-site for data acquisition. In addition, some electrical disturbances are rare in terms of location and frequency. Their measurement can therefore further prolong the data acquisition stage. Power quality monitoring devices can be placed permanently to cope with this issue. Nonetheless, processing and analyzing large amounts of data can also be time-consuming. Because of all these reasons, the automatic analysis of power quality data is a subject of growing demand, particularly for industrial and tertiary customers, due to the financial implications.

Problem statement and Objectives

This research work aims to develop an intelligent system to analyze power quality disturbances by processing electrical measurements (three-phase voltage and current) from a single monitoring device placed at the main electric supply point of an industrial or tertiary site.

Among the power quality disturbances affecting industrial grids, voltage sags are the most frequent and the most severe ones [START_REF]IEEE Recommended Practice for Monitoring Electric Power Quality[END_REF][START_REF] Wagner | Power quality and factory automation[END_REF][START_REF] Sarmiento | A voltage sag study in an industry with adjustable speed drives[END_REF], as they can cause equipment malfunction and unwanted stop of industrial processes. Therefore, they will be the topic of this research. The analysis of voltage sags is carried out in two steps: classification of voltage sag sources and analysis of their impact.

Understanding the origins of voltage sags is essential in the diagnosis process. The literature already proposes different solutions that report high accuracy for the classification of voltage sag causes. Our goal is to propose a solution that will have similar classification performance but will overcome several of their main limitations:

• The large majority of the state-of-the-art methodologies use statistical classifiers, which are highly dependent on data. However, access to power quality disturbance data, in the context of industrial grids, is not trivial.

• The generalization capabilities of these algorithms have not been evaluated in most of the cases. The accuracy of these methods mostly relies on training and testing data collected from the same source, which is not always practical in industrial applications.

• The physical and electrical interpretation of the fault features is weak. This makes the error analysis and trouble-shooting stage difficult to perform. Interpretable algorithms are preferable in industry for acceptability and reliability issues.

• There is a risk of information loss if only scalar features are extracted since the time dependence of electrical waveforms is not taken into account.

It should be noted that the algorithm we aim to develop has commercial purposes, hence its implementation must be compatible with industrial needs. This system should be able to be trained in factory using a reduced amount of synthetic data. The system should be able to perform effectively in different industrial sites without additional training. The decision-making process should be interpretable and it should provide a confidence index associated with the prediction to increase its acceptability and reliability.

Finally, once the source of a voltage sag is identified, its impact on the site should be assessed. We consider that a voltage sag has an impact if one or more loads trip after the voltage sag occurrence. The objective is to estimate the selfdisconnected load composition following a voltage sag. The goal is to take advantage of the dynamic nature of voltage sags to estimate the load composition. This study has its own challenges, but the most critical ones are:

• The access to a single monitoring location at the main electric supply point of the site.

• The lack of information on the industrial equipment downstream the monitoring point. The position of the monitoring device is fixed and imposed. The objective of the diagnosis system is to be as little intrusive as possible, with collected data from a unique monitoring point. This constraint limits the global cost of the solution with less monitoring material and reduces the risk of process dysfunction due to the installation and presence of external equipment on the industrial site.

• The diversity of electrical devices with different levels of sensitivity to voltage sags. The tolerance to voltage sags depends not only on the sag characteristics, but also on the characteristics of the industrial equipment. The development of a system capable of adapting to this equipment diversity is a real challenge.

To the best of our knowledge, the estimation of self-disconnected load composition after a voltage sag, with such constraints, has not yet been addressed in the literature. Providing as much information as possible to the industrial manufacturer about the characteristics of the disconnected devices after a voltage drop is essential for implementing adapted and cost-effective mitigation solutions, as only sensitive equipment should be targeted by these countermeasures.

Thesis outline

This PhD thesis is organized in five chapters:

• Chapter 1 introduces the basics and definitions of power quality. It also discusses the characteristics of industrial and distribution electrical grids in France. Additionally, it describes the main causes and consequences of power quality disturbances affecting industrial networks.

• Chapter 2 presents the state-of-the-art methods used in the analysis of power quality disturbances. Emphasis is given to the methods targeting the classification of voltage sags causes. A detailed comparison of the methods in the literature is presented, as well as a discussion on their main limitations.

• Chapter 3 presents a new method for the classification of voltage sag causes, given the voltage and current waveforms measured at the main energy supply point of the industrial site. The methodology follows a general fourstage scheme: data acquisition, pre-processing, feature extraction, and feature analysis. Each stage is presented in detail.

• Chapter 4 investigates the performance of the methodology presented in the previous chapter in terms of class separability, sensitivity to noise, sensitivity to fundamental frequency variations, and computational cost. We also evaluate the algorithm's global accuracy and generalization capabilities.

• Chapter 5 presents the problem of self-disconnected load composition estimation following a voltage sag. It introduces some of the most common methods in the literature in the domain of load estimation. Then, we propose an approach and apply it to a simple case study. We analyze the results obtained and based on these results we discuss the challenges and perspectives for future work.

Finally, a general conclusion and perspectives are provided.

Main contributions

The main contributions of this research work concern the development of a ready-to-implement voltage sag cause identification algorithm, based on the classification of multivariate time series signatures. Its main advantages are:

• The reduced amount of data necessary to build the reference signature database, which can be entirely composed of synthetic data.

• The provision of a confidence index associated with the prediction.

• The electrical intepretability of the signatures and the decision-making process.

• The robustness to noise levels up to SNR = 15 dB and to fundamental frequency variations up to ±0.5Hz.

• The good generalization capabilities when implemented on real field data, even for different industrial sites.

Finally, we also propose a first approach for the estimation of self-disconnected load composition. The preliminary results allowed to provide several guidelines for the improvement of the proposal in the future.

Chapter 1 Power Quality of Industrial Grids

Power quality analysis has drawn attention in the last few years due to the increasing sensitivity of electrical equipment integrating power electronics. A reliable energy supply ensures normal operation of the electrical equipment connected to the network. However, power quality disturbances can be responsible for disrupted production lines or services, equipment malfunction or even equipment damage, resulting into significant financial losses. Thus, the identification of the root causes of power quality disturbances and measure their effects on industrial equipment is a key step to provide adapted and cost-effective mitigation solutions, in order to reduce the impact on the site's productivity.

This chapter aims at providing relevant information on industrial electrical grids and power quality disturbances. Section 1.1 briefly describes some of the main characteristics of distribution and industrial grids, including the most representative categories of industrial loads. Section 1.2 presents the different power quality disturbances, as well as their impact in terms of financial losses for industrial customers. This section also explains the reason why this research work focuses on voltage sags, and the need to consider harmonic distortion in the analysis. Section 1.3 presents in detail the characterization of voltage sags, their main causes and consequences on industrial equipment, and Section 1.4 briefly introduces relevant notions regarding harmonic distortion, its main causes and consequences as well. Finally Section 1.5 closes the chapter.

Characteristics of Distribution and Industrial Grids

Distribution Grids in France

Distribution grids are intended to supply electricity to end consumers. In France, the standard voltage value for distribution networks is 20kV (medium voltage or MV). Only few portions of the grid still remain at 10, 15 and 33kV. The rated frequency is 50Hz, imposed by the transmission grid level.

Distribution substations are fed by the transmission network (63kV to 227kV) through HV/MV power transformers. All distribution feeders and substations are operated radially. Their rated power varies from 10 to 70 MVA. The total length of MV lines in France is about 622 187 km, consisting of overhead lines and underground cables. The MV network feeds more than 780 000 MV/LV substations. The presence of a grounding system at this level helps to reduce overvoltages when short-circuits occur by creating a low impedance path to earth. In France, two main grounding systems exist at the HV/MV substations : impedance and compensated neutral groundings. Fig. 1.1 illustrates a diagram of the distribution network.

Faults on line segments are common events in distribution and transmission networks. High currents due to short-circuits in the lines can seriously damage the power system equipment if they are not extinguished quickly. Thus, dedicated protection systems exist to protect the grid. Specific protection schemes are in place to minimize the impact on the power quality delivered to end customers. They relay on circuit breakers configured in a selective way to reduce the size of the impacted area, and automates (shunt and automatic circuit recolser) to benefit from faults' self-extinction.

High-consumption customers such as large industrial sites can be fed directly by the HV transmission network through a dedicated HV/MV transformer. Smaller industrial sites are fed by MV networks through MV/LV transformers and operate at low voltage levels (LV), set at 400 V/230 V (phase-to-phase and phase-to-ground voltage respectively).

LV industrial clients are supplied through MV feeders by a MV/LV transformer located on site, with rated power ranging from 50 to 1600kVA. The most common type of winding connections for transformers at this level is Dy. The wye-connection at the secondary allows an appropriate grounding system for the industrial facility, to reduce possible equipment damage due to faults. The most common grounding systems for industrial networks in France are the TN-C and

Main Industrial Loads and Equipment

There is a large variety of industries, types of processes and equipment. Depending on the activity sector, processes and loads can be very different. We can classify them in five categories. [START_REF] Repérant | Réseaux électriques industriels -Introduction[END_REF][START_REF] Mcgranaghan | Voltage sags in industrial systems[END_REF]:

1. Electrical motors. They are integrated in numerous types of processes as they convert electrical energy into mechanical energy, in order to power equipment. They represent more than 67% of the total electrical consumption in French industries [START_REF]La consommation d'énergie dans l'industrie en France[END_REF]. This category includes asynchronous motors, synchronous motors and DC motors. In particular, induction motors are the most widely used machines in industry due to their reliability and reduced cost. They include applications such as pumps, fans, compressors, conveyors, extruders, etc. 2. Power electronics. They are integrated in many industrial applications for improved control and energy efficiency. Their implementation has increased significantly in recent years. They include adjustable speed drives, AC/DC converters, soft-starters, etc. The most representative device of this category is the adjustable speed drive, which is coupled with induction motors for more flexibility in the speed control. 3. Electrothermal devices. They represent the heat-related equipment used in industry such as electrical heaters, induction/arc/resistance furnaces, welding machines, etc.

4. Electronic devices. They are used for control and surveillance of digital systems. They include programmable logic controllers (PLC), computers (PC), micro-controllers, etc. 5. Lighting. It includes different types of light-emitting devices such as classic incandescence lamps, gas-discharge lamps and light emitting diodes (LED).

Power Quality disturbances

The term power quality defined in [START_REF]IEEE Recommended Practice for Monitoring Electric Power Quality[END_REF] refers to a wide variety of electromagnetic phenomena that characterizes the voltage and current at a given time and at a given location in the power system. Power quality disturbances can also be defined as any change (distortion or fluctuation) in the power supply parameters outside the normative limits. These electrical disturbances can be caused by events in the transmission/distribution network (line faults, transformer and capacitor bank switching, etc.), or on the industrial site itself (large motor startup, power electronics, arc furnaces, etc).

Table 1.1 displays some of the electrical disturbances' characteristics according to their spectral content, duration and magnitude, as well as some of their main causes and consequences.

Impact of Electrical Disturbances on Industrial Grids

One of the main reasons for the rising research interest in the domain of power quality is the increasing implementation of power electronics in industrial processes due to their numerous benefits (more control flexibility, improved efficiency, reduction of energy losses, etc.). However, these electronic devices are sensitive to electromagnetic disturbances. At the same time, they generate disturbances, possibly affecting other loads. The term electromagnetic compatibility (EMC) [START_REF]IEEE Recommended Practice for Monitoring Electric Power Quality[END_REF] is used to describe the ability of a device to operate properly in accordance with the standards. In this sense, power quality can be seen as a means of ensuring EMC between the network and the connected loads.

Industrial processes are complex because they consist of different equipment such as motors, speed drives, contactors, digital control systems, etc. If a production unit contains equipment sensitive to a particular electrical disturbance, the whole process is at risk if such a disturbance occurs. The financial costs due to poor power quality can be particularly high, since the productivity and competitiveness in many industrial sectors depend on the continuity of the production. In [START_REF] Targosz | Pan-European power quality survey[END_REF], the economic losses in the industrial and service sectors of eight European 1 Frequencies in the range between 2-150 kHz are considered in the domain of supraharmonics [START_REF] Bollen | Standards for supraharmonics (2 to 150 kHz)[END_REF] countries were estimated up to €150bn, representing around 1% of their annual turnover, or more than 5% of their net profit.

Power quality costs can be divided into the following categories [START_REF] Chapman | The cost of poor power quality[END_REF]:

• Production interruption costs: some products or services are not available as long as the process is stopped, causing delay costs affecting the production. The duration of this interruption can vary from few seconds to few hours.

• Process restart costs: some industrial processes require several hours to be restored after an interruption. Several interventions may be necessary (cleaning, emptying, repairing, etc.) The costs include delay costs, human costs and material extra resources required to restart the process.

• Equipment damage costs: equipment can be damaged or can prematurely age following abrupt stop or overheating due to disturbances. Depending on the severity of the damage, reparation or complete replacement of several devices may be necessary, leading to additional costs.

• Defective products costs: some processes will not completely stop if a disturbance occurs. However, the high precision required for the fabrication of certain items can be affected resulting in defective products and leading to losses. The associated costs will depend on the possibility to recycle or discard such products.

• Other indirect costs, include penalties or compensations due to defective products or services. Mis-operation costs and energy losses due to inefficient process functioning can also generate financial losses to the customer. These costs are difficult to quantify. Depending on the industrial sector they can represent a significant percentage of the poor quality costs.

Among the power quality disturbances that mostly affect the most industrial power grids, are short interruptions (outages), voltage sags and harmonic distortion. Short interruptions along with voltage sags are the most severe disturbances with short-term consequences, since they can cause production interruption, defective products, process restarts, equipment damage and other indirect financial losses [START_REF] Targosz | Pan-European power quality survey[END_REF][START_REF] Wagner | Power quality and factory automation[END_REF]. Even if voltage interruptions have the highest impact, voltage sags are by far more frequent [START_REF]IEEE Recommended Practice for Monitoring Electric Power Quality[END_REF][START_REF] Wagner | Power quality and factory automation[END_REF][START_REF] Sarmiento | A voltage sag study in an industry with adjustable speed drives[END_REF]. Depending on the characteristics of the voltage sag and the sensitivity of the industrial devices, voltage sags can have different levels of impact on the customer. This is the main reason why this work focuses on the analysis of voltage sags. On the other hand, harmonic distortion is a steady-state disturbance with more long-term consequences. Due to the large presence of power electronic devices in industries, it is necessary to understand and consider the influence of harmonics in the developed analysis methods.

Voltage Sags

Characterization of voltage sags

A voltage sag is defined as a decrease between 10% and 90% in the RMS nominal voltage, with a duration between 0.5 cycle and 1 min [START_REF]IEEE Recommended Practice for Monitoring Electric Power Quality[END_REF]. The general agreement when referring to the amplitude of a sag is to take the remaining voltage as a percentage of the nominal voltage (IEEE standards 1346-1998, 493-2007, 1159-2019).

Three-phase voltage sags are usually characterized by their amplitude and duration. The amplitude of the sag is defined by the minimum voltage value within the three phases as illustrated in Fig. 1.2a. However this representation does not reveal all the characteristics of the threephase voltages. Since voltage (and current) are complex variables, each phase voltage is defined by a magnitude and a phase shift. The three phases can be represented in a phasor plot as shown in Fig. 1.2b. In a perfectly balanced threephase system, all voltage phasors have an amplitude of 1 pu and a phase shift of 120°between one another. During a voltage sag, the amplitudes and phases of the voltage phasors are affected (amplitude variation and phase-angle jump). These variations will depend on the characteristics of the sag, as it will be described in Section 1.3.2.

The time of occurence of the sag can also have an impact on some equipment. It is referred to as point-on-wave on sag initiation, and corresponds to the phase angle of voltage wave at which the voltage sag starts, taking as reference the zero crossing. Fig. 1.2c illustrates the point-on-wave (POW) of each phase at sag initiation. Although the impact of the POW is well-known for devices such as AC contactors, its effect is lower than the amplitude, duration or phase-angle jump. Voltage sags can be primarily classified as balanced or unbalanced, depending on the affected phases. However, there are different methods for characterizing and classifying voltage sags more precisely. One of the most widely used methods is the ABC classification proposed by Bollen et al. [START_REF] Bollen | Understanding Power Quality Problems[END_REF]. Voltage sags are described as a three phasor system, where the amplitude variations and phase-angle jumps for the three voltage components are represented. The method defines seven categories, each labeled from A to G. The ABC classification is presented by the author as a special case of the symmetrical-component classification method previously proposed in [START_REF] Zhan | Characteristic of voltage dips (sags) in power systems[END_REF].

Another popular method for the characterization of voltage sags is proposed by Ignatova et al. [START_REF] Ignatova | Space Vector Method for Voltage Dips and Swells Analysis[END_REF]. This technique is based on a space vector representation in the complex plane and the zero-sequence voltage for the classification of voltage dips and swells, using the Clarke transformation. Voltage sags are represented as ellipses in a 2d-space. Other techniques for the characterization of voltage sags can be found in [START_REF] Oubrahim | Disturbances Classification Based on a Model Order Selection Method for Power Quality Monitoring[END_REF][START_REF] Alam | Characterizing Voltage Sags and Swells Using Three-Phase Voltage Ellipse Parameters[END_REF][START_REF] Camarillo-Penaranda | Fault Classification and Voltage Sag Parameter Computation Using Voltage Ellipses[END_REF].

Main causes of voltage sags

The leading causes of voltage sags are: line faults, transformer energizing and large induction motor startups2 . These events significantly increase the absorbed current, creating a drop in voltage RMS values. The current and voltage waveforms depend on the event causing the voltage sag and on the loads connected to the grid. Distribution grid characteristics, particularly the available short-circuit power at the point of common coupling (PCC), has a significant impact on the voltage sag magnitude. These events can occur at the level of the transmission/distribution network, or on-site.

Voltage sags due to line faults

Line faults are the most common cause of severe voltage sags. They can occur in the power supply network or at the level of the industrial grid as well. They are caused by lightning strikes, wind, animals, or other objects in contact with energized lines (power supply network), or by short-circuits due to isolation failures or overloads (industrial grid).

Line faults can affect one or more phases. There are five fault types depending on the affected phases:

• Single phase-to-ground fault (LG) • Double phase fault (LL)

• Double phase-to-ground fault (LLG) • Three phase fault (LLL) • Three phase-to-ground fault (LLLG) They can also be classified in two groups: balanced and imbalanced faults. Balanced faults correspond to three-phase faults or three phase-to-ground faults, while imbalanced faults correspond to single-to-ground faults, double phase faults, and double phase-to-ground faults. As an example, a voltage sag caused by a double phase-to-ground fault (phases A and B) in the supply network is presented in Fig. 1.3. During the event, the faulted feeder causes a high inrush current, resulting in a voltage drop on the remaining feeders of the same bus. The voltage is restored as soon as the faulted feeder is detected and de-energized by the corresponding protective device (fault in the supply network), or as soon as the protection device or fuse trips due to overcurrent and disconnects the faulted asset (fault in the industrial site). Therefore, the duration of the sag depends on the settings and characteristics of the protection devices. For instance, in the case of faults in the power supply network, the close-open duty cycle of automatic circuit reclosers can cause two consecutive sags if the fault is not cleared successfully at the first reclose attempt.

The amplitude and phase-angle jump of the voltage phases during the sag depend on factors such as: the fault type, the fault resistance, the distance to the fault, the number and the type of transformers through which the sag propagates [START_REF] Bollen | Understanding Power Quality Problems[END_REF]. In particular, the propagation characteristics of the voltage sag in terms of amplitude and phase-shift will depend on the winding connections and the grounding system of the transformers between the fault location and the location of the measurement point [START_REF] Aung | The Influence of Transformer Winding Connections on the Propagation of Voltage Sags[END_REF].

Finally, it is important to point out the influence of directly connected induction motors since they contribute to the fault during the sag, as they operate as generators for a short period of time and reduce the voltage drop. This explains the non-rectangular shape of voltage sags in the RMS plot at the beginning of the sag [START_REF] Bollen | Understanding Power Quality Problems[END_REF][START_REF] Milanovic | The Influence of Induction Motors on Voltage Sag Propagation-Part I: Accounting for the Change in Sag Characteristics[END_REF]. Induction motors also have an impact on the voltage recovery, but this aspect will be detailed in Section 1.3.3.

Voltage sags due to a transformer energizing

In the distribution network, transformer energizing usually follows the end of a protection cycle. As the faulted feeder is disconnected, loaded transformers are energized after fault clearance. Transformer energizing can also take place in the industrial network if transformers are present at this level. This is particularly the case for MV sites, or LV sites with isolation transformers (used for protection of highly sensitive equipment).

During the transformer energizing, the magnetic flux in the windings can exceed the saturation limits, creating a high inrush current, thus a voltage sag. Transient magnetizing inrush currents can reach magnitudes as high as six to eight times the rated current [START_REF] Kulkarni | Transformer engineering: design and practice, ser. Power engineering[END_REF]. Core saturation affects the frequency spectrum of the waveforms during the event, introducing even harmonics, which are absent in normal operation. A significant level of voltage unbalance between the three phases takes place too. Another particular characteristic is the exponential voltage recovery, which is determined by the time constant related to the flux offset decay in the transformer core [START_REF] Bollen | Signal processing of power quality disturbances[END_REF]. Different parameters influence the characteristics of the sag caused by transformer energizing, although the most important are: the switching-on-angle at the energizing instant and the residual flux density in the transformer windings. Other factors such as the short-circuit power of the supply network and the load characteristics connected to the transformer influence the inrush current and the voltage recovery [START_REF] Kulkarni | Transformer engineering: design and practice, ser. Power engineering[END_REF]. Induction motor direct on-line startup induces a high inrush current that can reach 7 to 10 times the nominal value, resulting in a voltage sag at the level of the bus feeding the motor. As the voltage drops, the current suddenly increases. Then, the voltage recovers at the same rate, as the current slowly decreases. Fig. 1.5 illustrates the voltage waveforms and RMS values of a motor startup.

The total startup time depends on the motor electrical and mechanical characteristics, such as: rated power, inertia, electro-mechanical and load torque. It ranges from less than a second for small motors (5 to 45 kW), a few seconds for medium size motors (45 to 150 kW) and it can reach a few minutes for large motors (150 to 750 kW). We note that the typical power range for motors in LV industrial sites is less than 150 kW. An excessive sag may even prevent the motor from starting successfully [START_REF] Dugan | Electrical Power Systems Quality[END_REF]. This can be a serious issue when multiple motors are switched on at the same time.

Since the motor size and power requirements are known, the location and installation are planned in such a way that its startup does not affect other sensitive equipment. To limit the inrush currents, different starting methods are often used: autotransformer starters, resistance starters, delta-wye starters, soft-starters and variable speed drivers. Thus, voltage sags caused by induction motors are rarely deeper than 85%.

Effects and consequences of voltage sags

The impact level of voltage sags depends on the sag's characteristics, the grid conditions, and the connected equipment. Depending on the severity of the sag and the voltage tolerance characteristics of the processes, some will be able to ride through and continue normal operation. However, if the sag is severe enough to affect one or more sensitive devices involved in a certain process, the whole process can be interrupted.

Sensitive loads at risk of disconnection regarding voltage sags include: adjustable speed drives, electronic devices, induction motors, AC contactors and lighting [START_REF] Mcgranaghan | Voltage sags in industrial systems[END_REF][START_REF] Bollen | Understanding Power Quality Problems[END_REF]. The tolerance of a given device regarding the amplitude and the duration of voltage sags can be described in the so-called "voltage-tolerance curves" [START_REF] Bollen | Understanding Power Quality Problems[END_REF].

Voltage-tolerance curves

Electrical equipment is designed to operate under specific voltage conditions. If they are not respected, the devices are at risk of malfunction or simply may stop working. The voltage limit conditions of a device in terms of magnitude and duration are characterized in the voltage-tolerance curves (VTC). The first modern voltage-tolerance curve was introduced for computers by the Computer and Business Equipment Manufacturer's Association, known as the CBEMA curve [START_REF]CBEMA curve[END_REF]. A revised version and widely adopted reference known as the ITIC curve [START_REF]ITIC curve[END_REF], was published later by the Information Technology Industry Council (successor of CBEMA). The ITIC curve is intended to be used strictly on single-phase "information technology" equipment with 120V/60Hz-rated conditions, as testing in different conditions was not specified. More recently, the SEMI F47 [START_REF]SEMI F47 -Specification for Semiconductor Processing Equipment Voltage Sag Immunity[END_REF] curve has been proposed by the Semiconductor Equipment Materials International association for semiconductor processing equipment. Fig. 1.6 illustrates the tolerance limits defined by these three standard curves. Some noticeable differences can be observed in the area below 70% of voltage amplitude.

For other devices, the corresponding VTC is usually assimilated as a rectangular curve with a "knee" at V max and t max , as illustrated in Fig. 1.7. It defines two zones, a normal operation zone and a malfunction or trip zone. This curve helps to determine if a voltage sag with a given magnitude and duration will make the device trip or not. This, under the assumption that two voltage sags of the same magnitude and duration will have the same impact. However, this does not always hold true, especially for three-phase equipment, since three-phase voltage sags exhibit other characteristics than magnitude and duration, as described in Section 1.3.1. The influence of the other voltage sag attributes such as the type of sag, phase-angle jump and point-on-wave will depend on the nature of the affected device.

Induction motors

Directly fed induction motors have relatively good voltage-tolerance capabilities, although severe voltage sags can still stop them. The decay in voltage causes a drop of the electrical torque, which is proportional to the square of the voltage, thus a drop in the motor speed. The motor will be able to find a new operating point as long as the mechanical torque is higher than the load torque, otherwise it will continue to slow down until it stops. Induction motors usually can tolerate voltage sags of up to 70% of the remaining voltage [START_REF] Bollen | Understanding Power Quality Problems[END_REF].

Once the fault is cleared the voltage will recover. If the motor did not slow down completely, the voltage recovery will make the motor reaccelerate leading to an inrush current [START_REF] Gomez | Behavior of induction motor due to voltage sags and short interruptions[END_REF]. At this stage, motors draw a large current, up to 8 times their nominal current. Post-fault inrush current can lead to extended sags with a duration up to few seconds. If the dragged current is too high for too long, the overcurrent protection will activate and trip the motor, leading to a disconnection even after the actual voltage sag. On the other hand, if the motor dropped out during the sag but it is controlled by an ASD, an automatic reconnection operation can take place to restart the motor. Special cautions should be taken when reconnecting an induction motor that hasn't completely stopped. For instance, if the AC contactors are suddenly reclosed, there is a risk of "out-of-phase reconnection" between the supply voltage and the residual voltage in the motor circuits, leading to dangerously high inrush currents (up to 15 times nominal current) [START_REF] Baurand | La protection des moteurs BT[END_REF].

Adjustable Speed Drives

The configuration of three-phase adjustable speed drives (ASDs) consists of a rectifier, a DC bus, and a controlled inverter. When a voltage sag occurs, the rectifier provides a lower mean voltage level and the inverter is partly powered from the DC bus capacitor. If the sag is not too deep, an equilibrium point can be reached. However, if the voltage sag is lower than a threshold, the rectifier will stop conducting and it will only be able to supply the load for a short time, since the capacitor has a limited amount of energy. Normally, the undervoltage protection is set so it will activate before the minimum voltage threshold is reached [START_REF] Bollen | Understanding Power Quality Problems[END_REF].

In addition, ASDs' sensitivity highly depends on the type of voltage sag (balanced or unbalanced). Since the voltage delivered from the three-phase inverter to the DC bus depends on the phase-to-phase voltages, the ASD is able to ride through shallow unbalanced voltage sags. However, it is much more sensitive to balanced sags. It is also important to notice that the overcurrent transients due to the recharge of the DC bus at the end of the sag and voltage recovery, can be a potential source of damage [START_REF] Bollen | Understanding Power Quality Problems[END_REF]. To avoid this, an overcurrent protection is also often implemented. This means that the ASD can trip either during the voltage sag (due to the undervoltage protection, determining the duration limits), or at the voltage recovery (due to the overcurrent protection, determining the magnitude limits). Although ASDs voltage sag tolerance limits depend on various factors, they can generally tolerate sags of magnitude between 70-90% and of duration between 10-100 ms [START_REF] Djokic | Sensitivity of AC adjustable speed drives to voltage sags and short interruptions[END_REF].

Finally, phase-angle jumps caused by a sag can also be a source of malfunction since controlled inverters use phase-angle information (based on zero-crossings of the supply voltage waveform) for firing power switches.

Computers and PLCs

Computers, programmable logic controllers (PLCs), and other electronic devices are connected to the AC supply through a single-phase rectifier and a DC voltage controller. Similarly to ASDs, the tripping of these devices during a voltage sag depends on the DC bus capacity. Although their electrical consumption is much lower than other loads, the malfunction consequences due to voltage sags can be significant. Computers can suffer from information loss. PLCs, being a key component in complex processes, can produce incorrect control signals with serious damage to the process.

Since the inverters are single-phase, these devices are equally sensitive to bal-anced or unbalanced sags. Computers' voltage sag tolerance limits vary between 30-170 ms and 50-70%. Although extreme values can reach 8 ms, 88% and 210ms and 30% [START_REF] Bollen | Understanding Power Quality Problems[END_REF]. PLCs voltage tolerance characteristics are variable too, although they generally trip for sags of duration longer than 4 cycles (80 ms) and with a magnitude between 35-80%. These can be very sensitive equipment since some devices can trip as fast as after a half-cycle (10 ms) [START_REF] Bollen | Understanding Power Quality Problems[END_REF].

AC contactors

AC contactors connect motors to the power supply. Even if contactors are not high energy-consuming devices, they are key elements regarding voltage sags. The supply voltage powers an electromagnet which keeps the contact closed, and when the supply voltage drops under a certain threshold, the contact opens, preventing the motor from suddenly restarting when the voltage recovers. This mechanism protects the motor from unexpected and sometimes dangerous restarting. AC contactors usually trip for sags of magnitude lower than 30% and of duration superior to 10 ms, but in some cases they can trip at sags as shallow as 75%. The tripping of one motor caused by its contactor can turn down a whole process, even if the motor itself would have been able to ride through the sag.

In addition to voltage magnitude and duration, contactors are very sensitive to point-on-wave of sag initiation, as demonstrated by Djokic et al. in [START_REF] Djokic | Sensitivity of AC Coil Contactors to Voltage Sags, Short Interruptions, and Undervoltage Transients[END_REF]. The pointon-wave angle is measured between 0°and 90°, because the contactor's sensitivity exhibits a quarter-cycle symmetry. A voltage sag with point-on-wave at 90°will cause the tripping of the contactor more rapidly than another voltage sag of the same magnitude but with point-on-wave at 0°. This is because the current flowing through the coil is responsible for the electromagnetic force keeping the contactor closed. This characteristic explains the non-rectangular shape of the VTC curve for AC contactors. Finally, while voltage tolerance limits are mostly influenced by the POW [START_REF] Djokic | Sensitivity of AC Coil Contactors to Voltage Sags, Short Interruptions, and Undervoltage Transients[END_REF], they also depend on the phase-angle jump.

Gas-discharge based lighting

The voltage tolerance of lighting equipment depends on the type of lamp. In the case of gas-discharge lamps, voltage sags can extinguish the lamp, which needs to cool down for one to several minutes before restarting. They usually trip for sags of magnitude of 50% and 50ms duration, but these limits can also vary depending on their age. They become more sensitive to less severe sags as they need a larger voltage to operate. For instance, new lamps can tolerate sags of up to 45% magnitude, whereas older lamps at the end of their useful life will only tolerate sags of up to 85% magnitude [START_REF] Bollen | Understanding Power Quality Problems[END_REF].

Harmonic distortion

Harmonic distortion is due to the presence of nonlinear devices. Electronicbased loads are major contributors of harmonics in the power system. Due to the widespread of power electronics systems, harmonics have become a key issue in industrial, commercial and domestic installations [START_REF] Baggini | Handbook of Power Quality[END_REF].

Characterization of harmonics

Harmonics are sinusoidal voltages or currents whose frequencies are integer multiple of the fundamental frequency. The harmonic order n indicates the harmonic frequency f n , where f 1 =50Hz is the fundamental frequency. The distorted waveform corresponds to the superposition of one or more harmonic components to the fundamental frequency. Each harmonic component can be described by its amplitude, frequency and phase-shift. Fig. 1.8 illustrates the phasor nature of current harmonics, as well as the resulting waveform when they are superimposed.

Harmonic distortion levels in voltage or current waveforms at steady-state (periodic waveforms), can be characterized through the harmonic spectrum in the frequency domain, obtained with the Fourier Transform. This representation describes the amplitude and frequency of each harmonic. The amplitude of each component is given in percentage of the fundamental's amplitude.

Another way of characterizing harmonics is to use the Total Harmonic Distortion factor or T HD. The T HD factor for the current (T HD I ) is defined as in equation (1.1), where I 1 is the RMS value of the fundamental, I n represents the RMS values of the harmonic currents (f = n * 50Hz) up to the order N = 40 (harmonic currents over this limit are considered as in the supraharmonic domain). An equivalent formula is used for calculating the T HD U factor for voltage harmonic distortion as well.

T HD I (%) = N n=2 I 2 n I 1 * 100% (1.1)

Main causes of harmonic distortion

The main sources of harmonic static and transient currents in the industrial sector are [START_REF] Baggini | Handbook of Power Quality[END_REF][START_REF] Dugan | Electrical Power Systems Quality[END_REF]:

• Electronic devices, which nowadays are powered via switched-mode power supplies consisting of single-phase rectifiers. The current drawn by the rectifiers contains large amounts of harmonics. However their total power consumption is relatively low compared to other equipment (ie. induction motors).

• Three-phase converters, including variable speed drives and AC/DC con-verters mostly based on a three-phase diode bridge usually of "six-pulse" type, which produces current harmonics of order 6n ± 1 3,4 . Twelve-pulse bridges are also used in industry (although less common), and they mostly generate harmonics of order 12n ± 1.

• Fluorescent and LED lamps are highly non-linear. They are interfaced to the AC voltage supply through a single-phase rectifier. Their unitary power consumption is low, but since they are installed to lighten large surfaces, they can cause a noticeable level of harmonic distortion with mainly odd harmonics.

• Transformers do not induce a significant level of harmonics in normal operation. However, their harmonic content rises significantly within the saturation region. In particular during transformer energizing, transformers generate harmonics of order n = 2 and n = 3. The third harmonic is usually filtered by transformers of Dy type, which is a commonly used configuration in distribution and industrial networks.

• Arc furnaces produce harmonics of significant value from n = 2 up to n = 9. The current distortion is highly variable and changes depending on the heating stage of the furnace.

Table 1.3 illustrates the current waveforms and harmonic spectrum of some of the industrial devices responsible for introducing harmonic distortion into the industrial grid.

Effects and consequences of harmonic distortion

The injection of current harmonics by non linear loads to the grid may cause adverse effects on equipment. For instance the following devices are particularly sensitive to harmonic distortion [START_REF] Baggini | Handbook of Power Quality[END_REF][START_REF] Dugan | Electrical Power Systems Quality[END_REF]:

• Capacitor banks are often used by industrial facilities to improve their power factor in order to avoid financial penalties. However the presence of these devices can potentially magnify harmonic currents and cause resonance within the industrial grid. Resonance cause motors and transformers to overheat, and leads to the misoperation of sensitive electronic equipment. • Transformers experience additional heating, higher RMS current, eddy current and core losses. These additional losses may result in a higher operating temperature and premature aging.

• Motors normal operation can be affected by harmonic voltage distortion, which induces harmonic fluxes in the machine. These harmonic fluxes cause additional losses, decreased efficiency, heating and vibrations. These effects can be considered as mid-or long-term consequences, and if they are not corrected on time, they could reduce the motor's lifetime.

• Phase and neutral conductors, could overheat in the presence of current harmonics. Indeed, the total neutral current is almost zero in normal operation. However this is no longer the case in the presence of the third harmonic, which can add and flow through the neutral conductors. This can be problematic since neutral currents can reach high values, superior to the maximum capacity of neutral conductors.

• Converters and electronic devices can also be affected if they use controlled inverters based on zero-crossing information of the supply voltage waveform since harmonics can cause zero-crossing noise. Synchronization errors may occur and improper switching of semiconductor equipment can damage the device.

It is also important to note the effect of harmonic distortion in the calculation of power quantities, as described in Table 1.3. Under sinusoidal conditions (no harmonic distortion), the RMS voltage V and current I are equal to the RMS values of the fundamental frequency, V 1 and I 1 respectively. The apparent power S (VA), active power P (W) and reactive power Q (VAR) constitute a power triangle of angle cos(φ 1 ), equal to the phase-angle between voltage and current fundamental waveforms.

The calculation of the RMS current I and RMS voltage V taking into account harmonic distortion is given by the equation (1.2) and (1.3), respectively.

I = I 2 1 + N h=2 I 2 h (1.2) V = V 2 1 + N h=2 V 2 h (1.3)
Harmonic distortion mainly affects harmonic currents. Harmonic voltages are a consequence of harmonic currents flowing through the system. The voltage provided by the power supply will be distorted by harmonic currents depending on the magnitude of the source impedance. The smallest the source impedance, the lesser the impact of harmonic currents on the voltage. In most cases the source impedance is small enough to consider that the distortion of voltage waveforms is minimal. Thus, the voltage waveforms can be assumed to be sinusoidal.

In the presence of harmonic distortion, the power triangle does not hold anymore because the apparent power S contains cross term products due to the voltage and current harmonic components. Thus, the calculation of the apparent power S is modified by the introduction of a third power quantity, known as Distortion power or D (VAD), as described in Table 1.4. Finally it is also important to note that the power factor P F is equal to cos(φ 1 ) in the absence of harmonics, however this is no longer true in the presence of harmonic distortion.

Sinusoidal current Non-sinusoidal current

Active power

P = V 1 I 1 cos(φ 1 )
Reactive power

Q = V 1 I 1 sin(φ 1 )
Apparent power

S = V 1 I 1 S ≈ (V 1 I 1 ) 2 + V 2 1 N h=2 I 2 h Distortion power D = 0 D ≈ V 1 N h=2 I 2 h
Power factor P F = P S = cos(φ 1 ) P F = P S < cos(φ 1 )

Power triangle

S = P 2 + Q 2 S = P 2 + Q 2 + D 2

Conclusion

This chapter has introduced the electrical background necessary for the understanding and development of automatic analysis methods of power quality disturbances in industrial grids. Particular attention has been paied to voltage sags, which will be the focus of the research work developed in the following.

After the characterization of voltage sags in terms of amplitude, duration, phase-angle jump and point-on-wave, we have presented the main causes of voltage sags: line faults, transformer energizing and direct motor startup. In a similar way we have described the main effects and consequences of voltage sags on the most sensitive industrial equipment: induction motors, adjustable speed drives, computers and PLCs, AC contactors and gas-discharge lamps. The main causes and consequences of harmonic distortion were briefly introduced as well.

The following chapter presents the state of the art methods used for the analysis of power quality disturbances, including voltage sags.

Chapter 2

Power Quality Analysis: Application to Voltage Sags

Introduction

The purpose of this chapter is to present the methods and tools used in the literature for power quality (PQ) analysis. Section 2.2 describes the three main categories of PQ analysis: detection, location and classification1 of PQ disturbances, and Section 2.3 presents in more detail the automatic classification of PQ disturbances. A general five-stage scheme is presented, with a description of the state-of-the-art techniques. Then, in section 2.4 we focus on the classification of the causes of voltage sags. A summary of the most relevant references addressing this task is proposed, as well as a discussion on the main limitations of the methods proposed in the literature.

Analysis of power quality disturbances

Research on PQ analysis can be divided into three main categories: detection, location and classification of PQ disturbances, as illustrated in Fig. 2.1. Although the objective of each type of analysis is different, they share common methods. A brief description of the most relevant proposals for each category will be presented, with a focus on the classification of PQ disturbances.

Detection of Power Quality Disturbances

The first step in PQ analysis consists of detecting the presence of an anomaly in the voltage or current waveforms. It is usually based on a threshold analysis combined along with segmentation techniques to determine the disturbance starting and ending points. The precision to be achieved in this stage will depend on the requirements of the detection stage and the types of analysis to be performed (ie. location, classification, etc.). For instance, in [START_REF] Caujolle | Impact of Waveform Segmentation Accuracy on Disturbance Recognition Reliability[END_REF], the authors specified that a detection error higher than 1/4 of the duration of a cycle can have a significant negative impact on the calculation of features for the characterization of transformer energizing. Other authors aim for a disturbance detection precision within the first 1/16 cycle up to 1/8 cycle duration (or between 1.25 ms and 2.5 ms for a frequency of 50 Hz) [START_REF] Nagata | Real-time voltage sag detection and classification for power quality diagnostics[END_REF].

Waveform measurements can rapidly represent considerable amounts of data. Saving the totality of this data is not feasible due to memory storage capacity. For instance, a file containing the voltage and current waveforms during 1 s with a sampling frequency of 12.8 kHz, is about 1.5 MB in size, which corresponds to 5.4 GB for one single day. The transmission and storage of this data to computation centers can be financially costly as well. For this reason, PQ disturbance detection algorithms are intended to be implemented at the monitoring device level, in order to effectively detect the occurrence and only record the events of interest in detail. Since the detection stage is performed online, its computation time performance is essential. To cope with this, the methods and techniques for detection should be compatible with real-time computing.

Different strategies for disturbance detection have been proposed in the literature. For instance, Nagata et al. [START_REF] Nagata | Real-time voltage sag detection and classification for power quality diagnostics[END_REF] used Independent Component Analysis (ICA) for the detection and segmentation of voltage sags, obtaining near 100% detection accuracy for tests on data with SNR = 60 dB. But with results were significantly degraded with SNR = 40 dB. Bastos et al. [START_REF] Bastos | Universal Waveshape-Based Disturbance Detection in Power Quality Data Using Similarity Metrics[END_REF] proposed a general wave shape-based disturbance detection algorithm by comparing consecutive cycles of data using a similarity metric. The method is able to identify slight PQ disturbances without any assumptions on the system, and is robust to frequency variations. However, the detection window is one cycle long, and no precise segmentation within this one cycle window is performed.

Location of Power Quality Disturbances

Locating the origin of PQ disturbances (ie. voltage sags, transients, etc.) due to short-circuit faults, capacitor switching or transformer energizing is of great concern. The location methods can be divided into three types: relative location (upstream or downstream the monitoring point), area location (branch of the grid ) and exact location (precise distance from the monitoring point). The exact location of the source requires knowledge of the impedance and lengths of the electrical lines [START_REF] Stefanidou-Voziki | A review of fault location and classification methods in distribution grids[END_REF][START_REF] Pereira | Fault Location in Transmission Lines Using One-Terminal Postfault Voltage Data[END_REF][START_REF] Izykowski | Accurate Noniterative Fault Location Algorithm Utilizing Two-End Unsynchronized Measurements[END_REF]. Locating the source area requires knowledge of the topology of the network studied and the use of several measurement points [START_REF] Chang | An Efficient Approach for Tracking Transients Generated by Utility Shunt Capacitor Switching[END_REF][START_REF] Khani | A New Method for Online Determination of the Location of Switched Capacitor Banks in Distribution Systems[END_REF]. Relative location methods can be implemented using a single monitor point [START_REF] Parsons | A direction finder for power quality disturbances based upon disturbance power and energy[END_REF][START_REF] Shin | Signal Processing-Based Direction Finder for Transient Capacitor Switching Disturbances[END_REF].

Once the voltage dip has been detected in the industrial network, it must be established whether its origin is in the upstream network or the downstream network, i.e. at the industrial customer's site. This identification is essential to determine responsibilities and set financial penalties. We present hereafter a brief description from the literature of the methods for the relative location. Let us remind that one of the main constraints of our research work is the use of a single monitoring point and the absence of precise information on the industrial electrical installation and equipment. The methods for the relative location of the sources of voltage sags can be divided into five categories [START_REF] Mohammadi | Locating the source of voltage sags: Full review, introduction of generalized methods and numerical simulations[END_REF]:

1. Disturbance power and energy changes. Voltage sags are a consequence of events that can be considered as "energy sinks". The energy flow is analyzed to determine the source location. [START_REF] Parsons | A direction finder for power quality disturbances based upon disturbance power and energy[END_REF][START_REF] Kong | Voltage sag source location based on instantaneous energy detection[END_REF]. 2. Voltage-current characteristics. During a voltage sag, the current tends to increase for downstream events and to decrease for upstream events. The peaks or slopes of voltage and current trajectories are used to determine the location of the sag source [START_REF] Li | Method for voltage-sag-source detection by investigating slope of the system trajectory[END_REF][START_REF] Hamzah | A new approach to locate the voltage sag source using real current component[END_REF]. 3. Impedance changes. An estimated impedance is calculated, both in magnitude and angle. The variations of the sign and/or the angle are used to locate the direction of the source [START_REF] Tayjasanant | A Resistance Sign-Based Method for Voltage Sag Source Detection[END_REF][START_REF] Shao | Locating voltage sag source with impedance measurement[END_REF].

4. Voltage characteristics only. Using only voltage measurements, the voltage sag source direction can be estimated based on the phase angle jump and/or the change in magnitude [START_REF] Gomez | A novel methodology to locate originating points of voltage sags in electric power systems[END_REF]. 5. Current characteristics only. Sag source location is based on the changes of current phase and/or magnitude [START_REF] Pradhan | Fault Direction Estimation in Radial Distribution System Using Phase Change in Sequence Current[END_REF][START_REF] Polajžer | Instantaneous positivesequence current applied for detecting voltage sag sources[END_REF].

The majority of the methods proposed in the literature perform very well with balanced voltage sags in single-source radial networks. However, some of them tend to perform less effectively when applied to unbalanced sags or in meshed grids with distributed generation. Variants of these methods were proposed using positive-sequence phasors based on symmetrical components [START_REF] Paap | Symmetrical components in the time domain and their application to power network calculations[END_REF], providing notably improved results [START_REF] Mohammadi | Locating the source of voltage sags: Full review, introduction of generalized methods and numerical simulations[END_REF][START_REF] Polajzer | Evaluation of different methods for voltage sag source detection based on positive sequence components[END_REF]. Due to the transient nature of voltage sags, the instantaneous current positive-sequence [START_REF] Polajžer | Instantaneous positivesequence current applied for detecting voltage sag sources[END_REF] and the instantaneous current Clark components [START_REF] Polajžer | Detection of voltage sag sources based on the angle and norm changes in the instantaneous current vector written in Clarke's components[END_REF] proved to perform better compared to conventional phasor-based methods as presented in [START_REF] Mohammadi | Locating the source of voltage sags: Full review, introduction of generalized methods and numerical simulations[END_REF].

Notice that the effectiveness of the location methods depends on a proper detection of the event. If the beginning of the disturbance is not available or recorded, it will be impossible to determine the sag source location.

Classification of Power Quality Disturbances

In PQ analysis, there are three types of classification to be distinguished.

1. Classification based on the type of disturbance. The classes to be identified are: interruption, voltage sag, voltage swell, harmonic distortion, unbalance, flicker, etc. This domain has been widely addressed in the literature, and numerous classification methods have been proposed and reviewed [START_REF] Saini | Classification of power quality events -A review[END_REF][START_REF] Mishra | Power quality disturbance detection and classification using signal processing and soft computing techniques: A comprehensive review[END_REF][START_REF] Chawda | Comprehensive Review on Detection and Classification of Power Quality Disturbances in Utility Grid With Renewable Energy Penetration[END_REF][START_REF] Beniwal | A Critical Analysis of Methodologies for Detection and Classification of Power Quality Events in Smart Grid[END_REF]. Overall classification results are globally high. However, there are still some key issues to be addressed in this area. For instance, few proposals have been developed for three-phase systems and many methods are based and have been only tested on simulation data only. Although more recent methods have focused on the classification of single and multiple events, a generalized approach is still needed. Another aspect to be improved is the real-time performance which could be required for some applications. Finally, some authors have pointed out that more knowledge and electrical expertise should be incorporated into the algorithms, specially with statistical classifiers.

Classification based on the disturbance's type and characteristics.

Given a specific type of disturbance, the goal is to classify it according to its characteristics. For instance, voltage sags can be classified according to the affected phases, their magnitude and phase-angle jump. Several methods are based on characterization of voltage sags and swells [START_REF] Ignatova | Space Vector Method for Voltage Dips and Swells Analysis[END_REF][START_REF] Oubrahim | Disturbances Classification Based on a Model Order Selection Method for Power Quality Monitoring[END_REF][START_REF] Alam | Characterizing Voltage Sags and Swells Using Three-Phase Voltage Ellipse Parameters[END_REF], often using as reference the phasor-based ABC classification [START_REF] Bollen | Understanding Power Quality Problems[END_REF]. In this case, the classification's result can be observed directly from the measurements. Thus, these methods have the advantage of not requiring a training database for their development.

3. Classification based on the disturbance's cause. The objective of this type of classification is to determine the underlying event responsible of the observed disturbance. Indeed, every disturbance can be related to an event in the grid, but one event can induce different types of disturbances. For instance, a voltage sag can be caused by a line fault, a transformer energizing, or a starting motor. Steady-state harmonics can be caused by a variable speed drive, an arc furnace or an AC/DC converter. Transients can be caused by a transformer or capacitor switching.

Even though several techniques have been proposed in the literature based on the type and characteristics of the disturbance, few cause-based classification algorithms have been proposed. The need for methods that target the identification of the underlying causes of PQ disturbances has been pointed out by different authors [START_REF] Saxena | Power quality event classification: an overview and key issues[END_REF][START_REF] Khokhar | A comprehensive overview on signal processing and artificial intelligence techniques applications in classification of power quality disturbances[END_REF][START_REF] Mishra | Power quality disturbance detection and classification using signal processing and soft computing techniques: A comprehensive review[END_REF]. "A classification approach based on the underlying causes of disturbances may be more difficult but generally more relevant for the diagnosis of the electrical system", as highlighted by Bollen et al. in [START_REF] Bollen | Classification of Underlying Causes of Power Quality Disturbances: Deterministic versus Statistical Methods[END_REF]. Indeed, once a particular disturbance has been correctly detected and identified, it is essential to determine its cause and origin in order to apply appropriate countermeasures.

As presented in Chapter 1, this research work focuses on analysing voltage drops, in particular on the classification of their causes and determining their relative location. The following section presents the most relevant methods and techniques used in the literature for this type of classification, as well as their advantages and limitations.

Classification of power quality disturbances

The classical scheme of classification methods in the literature can be decomposed in the following five stages: data acquisition, pre-processing, feature extraction, feature selection and feature analysis, as illustrated in Fig. 2.2. In addition, it should be noticed that some tools can be used in more than one stage. For instance, Principal Component Analysis (PCA) can be used for feature extraction or for feature selection; meta-heuristic optimization techniques can be used in the feature selection stage or in the classification stage.

Data acquisition

This stage provides input data to the rest of the algorithm. Data can be either synthetic or real. Synthetic data can be obtained through mathematical equations or physics-based models describing the behavior of the studied system. Real measurement data on the other hand can be recorded experimentally either in a laboratory or in actual grids. However the main limitations of using measurements of actual grids, is the reduced number of disturbances that can be recorded, because they are operation anomalies.

Synthetic data generation

• Numerical models PQ disturbances can be synthesized using mathematical equations in order to mimic the behaviour of voltage or current waveforms. The simplest representation can be obtained using a pure sine single-phase voltage waveform, defined as in equation (2.1).

V (t) = Asin(wt -φ) (2.1)
Where A is the peak voltage amplitude, w = 2πf (f being the fundamental frequency) and φ the phase-angle. On this basis, we can model PQ disturbances in a simplified way as presented in Table 2.1 [START_REF] Igual | Integral mathematical model of power quality disturbances[END_REF]. Where u(x) is the step function, N the total number of cycles and T = 1/f . Numercial models are a more accessible source of data, especially for the early development of algorithms. They allow large amounts of data to be generated quickly while controlling the various parameters.

Table 2.1: Numerical modeling of power quality disturbances

PQ disturbance Mathematical model Sag V (t) = A[1 -α(u(t -t 1 ) -u(t -t 2 ))]sin(wt -φ) Swell V (t) = A[1 + β(u(t -t 1 ) -u(t -t 2 ))]sin(wt -φ) Interruption V (t) = A[1 -γ(u(t -t 1 ) -u(t -t 2 ))]sin(wt -φ) Harmonics V (t) = A[sin(wt -ϕ) + H h=2 α h sin(hwt -ϕ h )] Flicker V (t) = A[1 + λsin(2πf ′ t -ϕ)]sin(wt -φ) Transient impulse V (t) = A[sin(wt -φ) -ρ(e -750(t-ta) -e -344(t-ta) )(u(t -t a ) -u(t -t b ))] Parameters 0.1 ≤ α ≤ 0.9, T ≤ t 2 -t 1 ≤ (N -1)T 0.1 ≤ β ≤ 0.9, T ≤ t 2 -t 1 ≤ (N -1)T 0.9 ≤ γ ≤ 1, T ≤ t 2 -t 1 ≤ (N -1)T 2 ≤ H ≤ 40; ∀h ∈ [2, H], 0 ≤ α h 0.05 ≤ λ ≤ 0.1, 8Hz ≤ f ′ ≤ 25Hz 0.22 ≤ ρ ≤ 1.11, T ≤ t a ≤ (N -1)T , t b = t a + 1ms

• Simulation models

Although numerical modeling is a useful approach for data generation, it cannot reproduce all the dynamics that can be found in the disturbed waveforms. As presented in chapter 1, PQ disturbances are due to the response of power system components when specific events occur in the network (i.e. a voltage sag or swell caused by a line fault). In this sense, simulation models can generate more accurate waveforms as they emulate these same events and calculate the response of the power equipment based on the combination of more accurate physical equations. Some of the most widely used simulation software for the analysis of power systems are Matlab/Simulink, EMTP-RV and PowerFactory.

For instance, Fig. 2.3 illustrates the RMS2 voltages of a three-phase voltage sag obtained with the numerical equations versus a sag obtained with the simulation of a line fault close to an induction motor. The numerical model generates a perfectly square-shaped RMS voltage, as it cannot capture the influence of the induction motor during the voltage drop and recovery [START_REF] Bollen | Understanding Power Quality Problems[END_REF][START_REF] Milanovic | The Influence of Induction Motors on Voltage Sag Propagation-Part I: Accounting for the Change in Sag Characteristics[END_REF] as does the simulation model.

Real data measurement

Real data measurement is generally preferred over synthetic data as it naturally carries the physical phenomena and interactions, which otherwise would be difficult to reproduce with numerical or simulation models. Real data also allows evaluating the algorithms under real measurement conditions (different noise levels, frequency variations, etc). These measurements can be divided into experimental laboratory data and field data.

• Experimental laboratory data Data obtained in laboratory can be measured under controlled conditions. However, laboratory facilities can be oversimplified configurations depending on the availability of the equipment. For example, replicating an entire industrial site, even a small one, would require access to specific equipment and machinery. The generation of several types of disturbances simultaneously remains limited. Electrical interactions (resonance, harmonic distortion, etc.) between industrial loads and other network elements are more challenging to reproduce. This requires the use of strategies such as hardware-in-the-loop (HLP) simulation. This could be a good compromise between simulation and actual measurement data acquisition, provided that the models are sufficiently accurate.

• Field data

Field data are measured with PQ monitors placed at strategic locations in the studied network. The installation of such devices can be expensive depending on the hardware and software characteristics of the monitor. The installation can also require specific procedures such as power cuts depending on its location in the site's network, affecting the normal site's operation. To evaluate the performance of algorithms, only results obtained with field data are conclusive. However, this approach has several drawbacks, such as controlled measurement conditions and non-repeatability of events. The non-controlled conditions can induce incomplete information on the analyzed event. Similarly, the inherent nature of disturbances as anomalies can make it very difficult to generate a significant amount of data. This can negatively affect the constitution of a complete and meaningful data set, which is essential for developing and validating PQ analysis algorithms. This is not a problem for disturbance characterisation, since the ground truth can be extracted directly from the recorded waveforms. However missing information about the event conditions (source, connected devices, impact, etc.) can not always be extracted directly from the electrical waveforms.

Most of the methods in the literature use simulation data for their development and testing, as it is challenging to access large, representative and labelled real data sets. Simulation models also allow for generating data for a wide range of controlled parameters and for analysing the algorithm's performance in extreme scenarios that are rare and difficult to record in existing power systems. Nevertheless, validation with field data is mandatory if the algorithms are to be deployed in actual operating conditions. Therefore, a combination of simulation and field data seems to be a relevant approach to develop a robust solution.

Pre-processing

The pre-processing stage prepares raw data for further feature extraction and analysis. Pre-processing can consist of data cleaning, segmentation and format transformation.

Data cleaning

Errors can occur when measuring data. These errors are mainly due to the monitoring device, including missing data points, latched data, abnormal range or extreme values [START_REF]IEEE Recommended Practice for Monitoring Electric Power Quality[END_REF]. Tgese must be eliminated before further processing. Resampling data streams to obtain a more suitable sampling rate can be necessary. Denoising techniques can be applied to improve the overall performance of the algorithms. For instance, noise filters using wavelets such as Translation-Invariant Wavelet (TI-W) [START_REF] Liang | A translation-invariant wavelet representation algorithm with applications[END_REF] improved the classification accuracy in [START_REF] Mishra | Advanced signal processing and machine learning techniques for voltage sag causes detection in an electric power system[END_REF].

Segmentation

At the end of this stage, the electrical waveforms are segmented. Segmentation allows selecting only the section of the data stream to be processed. The required accuracy of the segmentation depends on the feature extraction methods. Segmentation methods with high precision are similar to those used for disturbance detection, presented in Section 2.2.1. When analyzing short-term disturbances, several feature extraction methods require the voltage waveforms to be split into transients and steady-state segments. A low accurate waveform segmentation can lead to a wrong calculation of features [START_REF] Caujolle | Impact of Waveform Segmentation Accuracy on Disturbance Recognition Reliability[END_REF], degradation of the global performance of the classification algorithm.

Format transformation

Electrical data are primarily measured as waveforms. This is the most usual format for electrical input data and it is generally preferred as it preserves the maximum information. Nonetheless, voltage and current waveforms can be transformed into other formats such as RMS values or images for further processing, as illustrated in Fig. 2.4. Normalization may be also applied at the end of this process.

Feature extraction

The obtained data after pre-processing can be used directly as input for the classification stage as some classifiers are designed to handle raw time series. Otherwise, a feature extraction stage is necessary. The goal of this stage is to apply one or more transformations and processing techniques to bring out the most relevant information from the input signal segments. A large majority of the methods in the literature apply one or more signal processing techniques and then extract numerical features from the processed signals. However, features are not restricted to numerical format. They can also be time series, or 2D and 3D images.

Time-frequency transformations

Time-frequency transformations are the most common methods used in the literature for PQ disturbance classification. They allow voltage (and current) waveforms to be represented in a time-frequency space.

• Short-time Fourier Transform Fourier transform (FT) [START_REF] Kehtarnavaz | CHAPTER 7 -Frequency Domain Processing[END_REF] is a well-known time-frequency transformation method, particularly useful for harmonic analysis. However, its implementation is limited to stationary signals. One of its variants, Short-Term Fourier Transform (STFT) [START_REF] Kehtarnavaz | CHAPTER 7 -Frequency Domain Processing[END_REF], can be applied to non-stationary signals. It provides information in the time-frequency domain as it computes the FT on short segments of the signal. STFT has been used for detection and feature extraction of PQ disturbances [START_REF] Bastos | Universal Waveshape-Based Disturbance Detection in Power Quality Data Using Similarity Metrics[END_REF][START_REF] Santoso | Characterization of distribution power quality events with Fourier and wavelet transforms[END_REF][START_REF] Youssef | Disturbance Classification Utilizing Dynamic Time Warping Classifier[END_REF]. STFT is also easy to implement. However, its main limitation is the time-frequency resolution, which depends on the length and overlap of the sliding time window.

• Wavelet Transform

Wavelet transform (WT) [START_REF] Heil | Continuous and Discrete Wavelet Transforms[END_REF] is one of the most popular time-frequency techniques for feature extraction. It has been applied for the analysis of PQ disturbances in [START_REF] Santoso | Power quality assessment via wavelet transform analysis[END_REF][START_REF] Chung | Power disturbance classifier using a rule-based method and wavelet packet-based hidden Markov model[END_REF][START_REF] Gaing | Wavelet-Based Neural Network for Power Disturbance Recognition and Classification[END_REF][START_REF] Silva | Fault Detection and Classification in Transmission Lines Based on Wavelet Transform and ANN[END_REF][START_REF] Erişti | Automatic classification of power quality events and disturbances using wavelet transform and support vector machines[END_REF], and more specifically for the classification of voltage sag causes in [START_REF] Wang | The method to reduce identification feature of different voltage sag disturbance source based on principal component analysis[END_REF][START_REF] Aggarwal | Recognition of Voltage Sag Causes using Vector Quantization based Orthogonal Wavelet[END_REF][START_REF] Saini | Fractionally delayed Legendre wavelet transform based detection and optimal features based classification of voltage sag causes[END_REF]. One of the advantages of WT over STFT is the possibility of having a variable resolution. Some of its variants are: Continous Wavelet Transform (CWT), Discrete Wavelet Transform (DWT) and Wavelet Packet Transform (WPT). DWT is usually coupled with a multiresolution analysis (MRA) and is less computational expensive than CWT. WPT can be interpreted as an extension of DWT, providing more precise frequency resolution than DWT, resulting in an equal-width sub-band filtering of the signals. WPT also has the advantage of being an orthogonal transform, which means that the energy in the signal is preserved.

Compared to the STFT, WT allows a fine temporal decomposition of the spectral content. But such a high level of decomposition can also be computational expensive. WT is also highly sensitive to noise. Finally, note that the performance of the transformation depends on the choice of the mother wavelet function. For PQ analysis, the daubechies 4 wavelet (Daub4) is usually implemented because of its good performance for disturbance detection and classification of short-term events, especially for slow frequency disturbances [START_REF] Brito | Daubechies wavelets in quality of electrical power[END_REF].

• Stockwell Transform

Stockwell transform (ST) [START_REF] Stockwell | Localization of the complex spectrum: the S transform[END_REF] is another time-frequency transformation tool, based on a scalable sliding Gaussian window. It can be interpreted as a generalization of STFT or an extension of CWT, overcoming some of their disadvantages as better resolutions in time and frequency are achieved. ST and its variants have been largely implemented for the classification of voltage sag causes [START_REF] Erişti | Automatic recognition system of underlying causes of power quality disturbances based on S-Transform and Extreme Learning Machine[END_REF][START_REF] Mishra | Advanced signal processing and machine learning techniques for voltage sag causes detection in an electric power system[END_REF][START_REF] Li | Adaptive S transform for feature extraction in voltage sags[END_REF], as well as for the detection and classification of other PQ disturbances [START_REF] Kumar | Recognition of Power-Quality Disturbances Using S-Transform-Based ANN Classifier and Rule-Based Decision Tree[END_REF][START_REF] Zhong | Power Quality Disturbance Recognition Based on Multiresolution S-Transform and Decision Tree[END_REF][START_REF] Tang | Classification of Complex Power Quality Disturbances Using Optimized S-Transform and Kernel SVM[END_REF]. However, besides the high computational cost, ST is not suitable for harmonic analysis because the widths of the frequency windows are directly related to their central frequency.

• Hilbert-Huang Transform

Hilbert-Huang transform (HHT) [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] is a method combining the Empirical Mode Decomposition (EMD) and the Hilbert transform (HT). EMD is a method used to recurrently decompose a signal into "modes" in the timedomain. The obtained intrinsic mode functions (IMF) can provide useful information contained in the signal, decomposing it into different spectral bands. Then, the HT is applied to the IMFs to obtain instantaneous frequency data, which makes it suitable for the analysis of non-stationary signals. Several authors have successfully implemented this technique in [START_REF] Afroni | Analysis of Nonstationary Power-Quality Waveforms Using Iterative Hilbert Huang Transform and SAX Algorithm[END_REF][START_REF] Biswal | Automatic Classification of Power Quality Events Using Balanced Neural Tree[END_REF][START_REF] Sahani | Automatic Power Quality Events Recognition Based on Hilbert Huang Transform and Weighted Bidirectional Extreme Learning Machine[END_REF] for disturbance type classification, and in [START_REF] Yang | Disturbance Source Identification of Voltage Sags Based On Hilbert-Huang Transform[END_REF][START_REF] Manjula | Empirical mode decomposition with Hilbert transform for classification of voltage sag causes using probabilistic neural network[END_REF] for voltage sag cause classification. One advantage of HHT is that the composition does not require a window selection. On the other hand, it has limited performances for discriminating components in narrow band signals, and the time-frequency decomposition is more difficult to interpret than other methods.

• Variational Mode Decomposition

Variational Mode Decomposition (VMD) [START_REF] Dragomiretskiy | Variational Mode Decomposition[END_REF] is a variant of EMD. It was developed in an attempt to improve the performance of EMD, in particular in terms of sensitivity to noise. Mishra et al. [START_REF] Mishra | Advanced signal processing and machine learning techniques for voltage sag causes detection in an electric power system[END_REF] used a combination of ST and VMD for the classification of voltage sag causes, and Sahani et al. [START_REF] Sahani | Variational mode decomposition and weighted online sequential extreme learning machine for power quality event patterns recognition[END_REF] used VMD for the identification of single and multiple PQ disturbances. Although VMD is more robust to noise, its main drawback is the boundary effects.

Space transformation methods

Space transformation methods are more rarely used for the classification of PQ disturbances according to their type or cause. Their implementation mainly concerns the characterization and localisation of disturbances.

• Clarke transform

Clarke transform is used to transform a three phase system (X A , X B , X C ) into an orthogonal system of three components (X α , X β , X γ ). In a balanced system, the third component X γ is zero, which simplifies the analysis. It is a power invariant transformation particularly used for space vector modulation control in AC drives.

In [START_REF] Ignatova | Space Vector Method for Voltage Dips and Swells Analysis[END_REF], the authors used a space vector method derived from the Clarke transformation using the ellipse-shaped Clarke's components when plotted in a cartesian system for voltage sag and swell characterization. A similar approach was used in [START_REF] Camarillo-Penaranda | Fault Classification and Voltage Sag Parameter Computation Using Voltage Ellipses[END_REF] for the characterization of fault-caused voltage sags. In [START_REF] Polajžer | Detection of voltage sag sources based on the angle and norm changes in the instantaneous current vector written in Clarke's components[END_REF], the Clarke's components of the instantaneous line current vector was used for the localisation of the sources of voltage sags. Clarke transform is a time-domain analysis tool, which means that no information in the frequency domain is provided.

• Fortescue transform

Fortescue transform [START_REF] Fortescue | Method of Symmetrical Co-Ordinates Applied to the Solution of Polyphase Networks[END_REF] is a linear transformation in the complex domain widely used in the analysis of unbalanced three-phase power systems. It allows representing an unbalanced set of three phasors (X A , X B , X C ) into a new balanced system known as symmetrical components (X + , X -, X 0 ).

In [START_REF] Bollen | Different methods for classification of three-phase unbalanced voltage dips due to faults[END_REF][START_REF] Oubrahim | Disturbances Classification Based on a Model Order Selection Method for Power Quality Monitoring[END_REF] the authors proposed an algorithm for voltage sag characterization based on the analysis of the symmetrical components. As mentioned in section 2.2.2, the use of positive-sequence phasors improved the results of voltage sag source location methods in [START_REF] Polajzer | Evaluation of different methods for voltage sag source detection based on positive sequence components[END_REF] when applied to asymmetrical sags, and in [START_REF] Polajžer | Instantaneous positivesequence current applied for detecting voltage sag sources[END_REF][START_REF] Mohammadi | Employing instantaneous positive sequence symmetrical components for voltage sag source relative location[END_REF] the authors went a step further by employing instantaneous symmetrical components for the same task. A symmetrical component-based modified technique has also been used in [START_REF] Kumar | Symmetrical Components-Based Modified Technique for Power-Quality Disturbances Detection and Classification[END_REF] for the classification of disturbance type. As for Clarke transform, Fortescue transform is a time domain technique with no spectral information provided. Nevertheless, the main advantage of Fortescue transform over other transformations is the electrical interpretability of its symmetrical components. Indeed, different electrical power equipment (generators, lines, transformers, motors, etc.) can be modeled using equivalent circuits based on symmetrical components [START_REF] Das | Understanding Symmetrical Components for Power System Modeling[END_REF]. Some of the main benefits and drawbacks of the presented signal processing techniques are synthesized in Table 2.2. However, a detailed comparison of their performance and efficiency is not provided since there is no standard open source PQ disturbance database for benchmarking. As a consequence, each author uses custom data sets, which differ in size, source (synthetic or real), noise levels, categories of PQ disturbances analyzed and classification goal (type-, characteristicor cause-based).

Signal descriptors and statistical moments

Statistical parameters and signal descriptors are numerical values. They are extracted from transformed signals and grouped into a feature vector, which is frequently used as unique output of the feature extraction stage.

• Signal descriptors

They include various numerical features that are calculated from time series signals such as: minimum and maximum values, energy, spectral entropy, Shannon entropy, etc. They also include electrical descriptors such as total harmonic distortion (THD).

• Statistical moments

Statistical moments are frequently used for characterizing the transformed signals and their subcomponents. They include mean, standard deviation, also called low-order statistics, and skewness and kurtosis, also referd as high-order statistics (HOS).

Feature selection

Feature selection is an optional stage between feature extraction and classification stages. The goal is to optimally choose the most relevant features for the final feature analysis stage. Indeed, feature extraction can result in a large number of features including redundant and correlated values. The complexity of the algorithms can rapidly increase if the number of features is too large. In fact, the performance of many classifiers such as machine learning models highly depend on the quality of the selected input features. Although the feature selection stage is not systematically included in the classification methods, some authors have

VMD

Relative energy ratio (RER), mode instantaneous amplitude (IA), number of zero crossings, center frequencies from each of the modes (four levels) [START_REF] Mishra | Advanced signal processing and machine learning techniques for voltage sag causes detection in an electric power system[END_REF].

pointed out the importance of selecting relevant and non redundant features and their impact on the classifiers' results [START_REF] Gunal | The search for optimal feature set in power quality event classification[END_REF].

Filters, wrappers and embedded methods

There are three classic categories for feature selection algorithms: filter-based, wrapper-based and embedded-based methods [START_REF] Chandrashekar | A survey on feature selection methods[END_REF].

• Filter methods define the relevance of features (or feature importance) using variable ranking techniques as criteria for selection, ie. Pearson correlation, Chi-square test, etc. The selection does not depend on the classification's performance. In [START_REF] Hafiz | Feature selection for power quality event identification[END_REF], a filter-based approach was used for PQ event identification, achieving a feature size reduction between 50 to 86% without any compromise in the classification performance. The main advantage of filters is the possibility to obtain a generic feature subset independently of the chosen classifier. For this reason they are faster than wrappers and embedded methods and it is a good approach when the number of features is large. It avoids overfitting but it may fail to effectively perform the best selection of features.

• Wrapper methods use the classification output as the objective function to evaluate the optimal feature subset. The subsets of features are obtained through search algorithms such as sequential selection methods. In [ [START_REF] Biswal | Detection and characterization of multiple power quality disturbances with a fast S-transform and decision tree based classifier[END_REF] incorporates the feature selection in the process of design of a DT classifier using Gini index as optimization criterion. Embedded methods are more computationally expensive than filters but they are faster than wrappers. They are also very effective and less prone to overfitting. However, they cannot be associated with any classifier, as the feature selection process needs to be integrated into the classifier structure.

Metaheuristic optimization techniques are also popular for feature selection. They are usually implemented as wrappers using the classifier output as objective function, but they can also be used in a filter-based approach if the objective function is defined with a metric independently of the classifier. Search algorithms include Genetic Algorithms (GA), Particle Swarm Optimization (PSO), Ant Lion Optimization (ALO), Teaching-learning-based Optimization (TLBO), etc. These nature-inspired algorithms can be applied for the optimal selection of features, but also for the classifiers' parameters optimization. Jamali et al. [START_REF] Jamali | Identification of optimal features for fast and accurate classification of power quality disturbances[END_REF] use a combination of GA, SFS and mRMR for optimal selection of various sets of features. Saini et al. [START_REF] Saini | Fractionally delayed Legendre wavelet transform based detection and optimal features based classification of voltage sag causes[END_REF] implement ALO for feature selection applied to voltage sag cause classification. Behera et al. [START_REF] Behera | Power quality time series data mining using S-transform and fuzzy expert system[END_REF] use a fuzzy adaptive PSO technique for optimizing expert system rules for PQ disturbance type classification.

Dimensionality reduction methods

• Principal Component Analysis

Principal component analysis (PCA) is a statistical technique for feature dimensionality reduction that allows the representation of a set of features into a new space of linearly uncorrelated features. Reducing the dimension of the feature vector allows to reduce the redundancy of features and to better control the complexity of the classification algorithms, also known as the "curse of dimensionality". In [START_REF] Wang | The method to reduce identification feature of different voltage sag disturbance source based on principal component analysis[END_REF] PCA is applied to reduce the number of features, improving the voltage sag source overall classification results when using a Support Vector Machine (SVM) as classifier. The use of a reduced subset of features decreased the computational cost and significantly improved noise sensitivity of the SVM (up to 40%) compared to the classifier trained with the original feature set.

• Autoencoder

Autoencoders are a type of unsupervised neural network, used for learning data representation while reducing dimensionality. The dimensionality reduction is learned and data-driven. They are composed of two subnetworks: an encoder and a decoder. The encoder compresses the input data into a latent-space representation, and the decoder is able to reconstruct the original input using the latent-space representation. The autoencoder is trained using both subnetworks, although the encoder is in reality the structure of interest as it creates the latent-space representation. This type of architecture is usually embedded inside deep neural networks. Autoencoders are used for unsupervised feature extraction and selection in [START_REF] Li | Self-Supervised Voltage Sag Source Identification Method Based on CNN[END_REF], as part of a deep neural network model for voltage sag source's classification.

Feature analysis

The final stage of the process is feature analysis. There are two main categories: rule-based and machine learning techniques. Rule-based methods include threshold analysis and expert systems. The most recent literature proposals use machine learning techniques for classification, as they are more flexible and achieve better overall results. However, the complexity level and the performance of these classifiers are usually correlated to the amount of necessary training data, which varies between 100 to 500 events per class in the reviewed literature. Some of the techniques in both categories have also been combined with other approaches such as Fuzzy Logic3 , giving as result hybrid methods such as fuzzy expert systems and neuro-fuzzy systems.

We present here a brief description of the most common methods in the literature to provide a global overview. Table 2.4 summarizes their main benefits and drawbacks. As for feature extraction methods in Section 2.3.3, a detailed comparison of their performance and efficiency is not provided since the data used for training and testing, as well as the number of classes and classification goals differ from one publication to another. However, such detailed results are presented in Section 2.4, where methods applied to the classification of voltage sag causes are addressed and can be effectively compared as in Table 2.5.

Rule-based methods

• Threshold analysis

One of the most basic techniques for decision-making is threshold analysis. The threshold values can be defined using theoretical values or calculated values from observations of the data. In [START_REF] Ding | Voltage Sag Disturbance Detection Based on RMS Voltage Method[END_REF] electrical features extracted from RMS voltages are used (amplitude, duration, phase-asymmetry, etc.), and the voltage sag source identification is performed based on set of rules according to pre-defined thresholds. Similarly in [START_REF] Nunez | Feature characterization of power quality events according to their underlying causes[END_REF], electrical features such as changes in active and reactive power, second order harmonic amplitude, RMS voltage shape, etc. are used to define a rule-based framework. In other cases, thresholds are defined for features obtained after processing the original waveforms (Hilbert-Huang transform [START_REF] Yang | Disturbance Source Identification of Voltage Sags Based On Hilbert-Huang Transform[END_REF], and Clarke transform [START_REF] Ignatova | Space Vector Method for Voltage Dips and Swells Analysis[END_REF]). When thresholds are defined through data observation, their robustness depends on the size and representativity of the analyzed data. The advantages of rule-based methods through threshold analysis is their easy implementation and that no learning process is needed. Depending on the method used for determining the thresholds, generalization capabilities can not be guaranteed.

• Expert Systems

An expert system is an algorithm that represents knowledge and mimics the decision-making ability of a human expert. A set of rules is defined based on this knowledge and its performance depends strongly on human expertise. It usually consists of a user interface, an inference engine and a knowledge base. In [START_REF] Styvaktakis | Expert system for classification and analysis of power system events[END_REF] an expert system is used for the classification of different types of power system events causing voltage sags and interruptions. Similarly, a fuzzy expert system is employed in [START_REF] Caujolle | Identification et caractérisation des perturbations affectant les réseaux électriques HTA[END_REF] for the identification of power system events in MV grids. As for threshold analysis, expert systems do not require large amounts of data since there is no learning process. Their modular structure allows to easily add new classes if needed, without requiring to modify the existing blocks (or re-training as it is the case for statistical classifiers). However, they are difficult to develop and require expert knowledge.

Machine learning methods

Machine learning methods are subdivided in different subcategories. They can be supervised, unsupervised or semi-supervised, depending on the characteristics of the input data. In supervised learning, a set of labeled data is given to the algorithm for "training" the model's internal parameters to better predict the output labels. In unsupervised learning the data is unlabeled, and in semi-supervised learning only a small amount of data is labeled, combined with a large amount of unlabeled data. Depending on the final task, machine learning models are also divided into classification or regression models. The goal of classification models is to predict a categorical label, whereas regression models are used to predict continuous numerical values.

In the reviewed literature, the totality of the machine learning methods used for the classification of PQ disturbances belong to the subcategory of supervised classification models.

• Decision tree

A Decision Tree (DT) is a machine learning method that has a flowchart-like structure. Each internal node represents a condition on a specific feature, each branch represents a possible outcome based on the condition and each leaf node represents a class label. The method is based on the choice of a feature that maximizes and fixes data division. There are different criteria used to select the feature that maximizes data division, such as Gini index, entropy, information gain, etc. In [START_REF] Biswal | Detection and characterization of multiple power quality disturbances with a fast S-transform and decision tree based classifier[END_REF] this technique is used in a solution proposed for the classification of PQ disturbances. The Gini diversity index is used as fitness measure for the construction of the decision tree. Compared to manually defined thresholds in rule-based techniques, DTs' decision rules are obtained automatically from features. They also have the advantage of having an interpretable decision-making process. The main drawback is their instability risk, as small changes in data can cause large structural changes.

• Random Forest Random Forests (RF) consists of many DTs combined together to get a more precise result. It is less prone to overfitting, and is generally more robust. However it is more computationally costly and more difficult to visualize and interpret compared to DTs. In [START_REF] Liu | Identification and Location of Voltage Sag Sources Based on Multi-Label Random Forest[END_REF][START_REF] Yikun | Classification Method of Voltage Sag Sources Based on Sequential Trajectory Feature Learning Algorithm[END_REF], the authors use a RF classifier for the classification of voltage sag causes. They obtained better classification results with RF when compared to SVM, ANN and DTs, using the same set of features (between 3% to 7% increase in global accuracy).

• Support Vector Machine

A Support Vector Machine (SVM) is a supervised model-based on statistical learning theory. Its objective is to construct an optimal decision function that accurately predicts new data by minimizing the classification error. The separation of classes is achieved by minimizing the margin of separation between classes in a high-dimensional feature space where the initial features have been implicitly mapped. Feature mapping into a new high-dimensional space is performed through kernel functions. There are different kernel functions: Linear, Polynomial, Gaussian, Radial Basis Function (RBF), Sigmoid, etc. Other parameters to be tuned include the regularisation parameter C, that allows to control the trade-off between misclassification tolerance and overfitting, and the parameter γ that regulates the influence of the most distant data points on the decision boundary definition. The choice of an adapted kernel function and correct parameter definition have a significant effect on the efficiency of the SVM.

SVM is a technique frequently used for PQ disturbance classification tasks. In [START_REF] Axelberg | Support Vector Machine for Classification of Voltage Disturbances[END_REF][START_REF] Wang | The method to reduce identification feature of different voltage sag disturbance source based on principal component analysis[END_REF], the authors implement a SVMs and in [START_REF] Sha | Identification Method for Voltage Sags Based on K-means-Singular Value Decomposition and Least Squares Support Vector Machine[END_REF] Sha et al. proposed a Least Squares SVM (LS-SVM) for the classification of voltage sag causes. One of its main benefits is the good trade-off between generalization performance and complexity. It is also fast to train compared to other approaches. However its performance depends highly on the kernel selection, and the final model is difficult to interpret.

• Artificial Neural Network

An Artificial Neural Network (ANN) is an statistical model capable of learning non-linear functions. An ANN is typically composed of an input layer, an output layer, and one or more hidden layers. ANN is the most most popular classification technique used in the classification of PQ disturbances. The most implemented variant is the multilayer perceptron (MLP). The inner neuron's parameters (weights and biases) are optimized in an iterative way using a gradient descent-based method, which is usually backpropagation.

MLP is implemented in [START_REF] Nagata | Real-time voltage sag detection and classification for power quality diagnostics[END_REF] for the classification of voltage sag sources. The authors report better classification results using MLP compared to SVM (+7%), where accurate results are obtained even with features calculated on the first 1/2 cycle window. ANNs are versatile and are able to represent complex non-linear functions. However they are highly data-driven, difficult to interpret and prone to overfitting. Other popular variants used for PQ analysis include Probabilistic Neural Network (PNN), Extreme Learning Machine (ELM) and Neuro-fuzzy Systems (NFS).

PNN is a neural network based on Bayesian theory. Unlike other neural networks, PNNs do not require a learning process or initial weights. Manjula et al. [START_REF] Manjula | Empirical mode decomposition with Hilbert transform for classification of voltage sag causes using probabilistic neural network[END_REF] reported better results with a PNN as classifier compared to MLP (+3%).

ELM [START_REF] Huang | Extreme learning machine: Theory and applications[END_REF] is based on a single hidden layer feedforward neural network. Different from gradient-based methods, ELM assigns random values to some of its weights which are frozen during training. This technique allows the algorithm to converge much faster than gradient-based ANNs, while achieving good generalization performance. ELM shows higher accuracy when compared to ANN, k-Nearest Neighbors (kNN) and SVM in [START_REF] Mishra | Advanced signal processing and machine learning techniques for voltage sag causes detection in an electric power system[END_REF] and better trade-off between accuracy, computational cost and noise sensitivity when compared to PNN and two different fuzzy clustering classifiers in [START_REF] Li | Adaptive S transform for feature extraction in voltage sags[END_REF].

Finally, NFS is a combination of ANN and fuzzy logic. Their fast learning and generalization capabilities are their main benefits. In [START_REF] Pires | Power quality disturbances classification using the 3-D space representation and PCA based neurofuzzy approach[END_REF] a 3D space and PCA-based approaches are combined with NFS for the classification of various PQ disturbances.

• Deep Neural Network

In recent years deep learning (DL) algorithms have gained attention from different fields such as speech recognition, computer vision and signal processing. DL algorithms have the ability to learn optimal features from training data, avoiding manual extraction and selection of features. Deep Neural Networks (DNN) contain multiple layers resulting in higher complexity during the training stage, but can outperform classical ANNs. Increasing the depth (number of layers) increases the capacity of the model. However if the model is too deep and the training dataset is not large enough, it can rapidly lead to overfitting and poor generalization capabilities.

Some of the most commonly used DNN architectures include Convolutional Neural Networks (CNN) and Long-Short Term Memory networks (LSTM).

CNNs use convolution for processing and extracting features of 2D or 3D data, one of the reasons why it is mainly used for image processing. In [START_REF] Li | Self-Supervised Voltage Sag Source Identification Method Based on CNN[END_REF][START_REF] Zheng | Recognition method of voltage sag causes based on two-dimensional transform and deep learning hybrid model[END_REF], the voltage waveforms are converted to 2D gray-scaled images and are fed to the network as input. LSTM is a type of recurrent neural network (RNN). The feedback connections allows the model to learn dependencies between sequence data, making them more suitable for the analysis of time series. In [START_REF] Zheng | Recognition method of voltage sag causes based on Bi-LSTM[END_REF][START_REF] Wang | Method of Voltage Sag Causes Based on Bidirectional LSTM and Attention Mechanism[END_REF] the authors successfully implement a bidirectional LSTM model (Bi-LSTM) for voltage sag source identification. Some of the main properties of the presented classification techniques are summarized in Table 2.4.

Classification and relative localisation of voltage sag causes

In the previous section, we presented a general scheme for the classification of PQ disturbances based on their type, characteristics or underlying causes. This scheme as well as the presented methods and techniques are also valid for the classification of voltage sag causes. Table 2.5 summarizes the main characteristics of the most relevant proposals in the literature. All the methods are capable of identifying the three main causes of voltage sags: line faults (F), transformer energizing (TE) and induction motor startup (MS). Some of these proposals are also able to identify combinations (C) of the previous categories (ie. voltage sag caused by a line fault followed by a motor starting). Other methods are designed to differentiate multiple types of short-circuit faults depending on the affected phases. They are noted as F(n), n being the number of fault types.

Discussion and main limitations of the methods in the literature

The global classification accuracy of the methods referenced in Table 2.5 is very high. Nevertheless, there are some limitations in the presented approaches regarding training and testing data that should be highlighted. The major ones can be considered as :

1. Large amount of labeled training data is required. Most of the approaches in the literature are based on statistical classifiers, thus highly data-driven. Machine learning techniques, and particularly deep learning networks, require a significant amount of data for effective training. Constituting a large and diverse dataset can be challenging since voltage sags are anomalies, and therefore not frequent events. In addition, the dataset must be labeled as classifiers are supervised, which can be difficult to obtain.

2. The generalization capabilities have not been fully evaluated. The large majority of the methods need to be trained and tested using data from the same data pool (either simulation or field data). Only [START_REF] Zheng | Recognition method of voltage sag causes based on two-dimensional transform and deep learning hybrid model[END_REF][START_REF] Wang | Method of Voltage Sag Causes Based on Bidirectional LSTM and Attention Mechanism[END_REF][START_REF] Zheng | Recognition method of voltage sag causes based on Bi-LSTM[END_REF] present results using a mix of simulation and field data. Although the classifier is firstly trained with synthetic data, it still requires a portion of field data for the final parameter tuning. The authors explain that this strategy is used with the objective of accelerating the learning process and improving overall results. However, results of these networks exclusively trained with simulation data and tested on field data have not been reported.

Nevertheless, we know from [START_REF] Bollen | Classification of Underlying Causes of Power Quality Disturbances: Deterministic versus Statistical Methods[END_REF][START_REF] Axelberg | Support Vector Machine for Classification of Voltage Disturbances[END_REF] that this is not a trivial task. Both authors report good results when training a SVM classifier with field data from one network and testing it on a different one. However, the classifier's performance significantly deteriorate when trained with simulation data and tested with real data. The same authors also point out that the usefulness of these methods for commercial purposes depends on the capacity of the classifier to be pre-trained on factory and work accurately when deployed in different networks, as it is not realistic that customers can train the algorithms themselves. In addition, from an implementation point of view, it is highly demanded that classifiers are purely based on synthetic data, which is easier to obtain [START_REF] Axelberg | Support Vector Machine for Classification of Voltage Disturbances[END_REF].

In a similar way, the performance of the algorithms depend on the choice of features, which should be done carefully to avoid the loss of relevant information. Thus, there are some aspects regarding the feature extraction process to be noted:

1. Low interpretability of the extracted features. Most of the used transformations do not integrate the physical and electrical properties of the PQ events. Thus, trouble-shooting and error analysis become tedious since most machine learning classifiers are "black-box" models. The generalization capabilities of the algorithms are also compromised, as it is difficult to provide guarantees on the behaviour of the selected features when applied to new field data. Therefore, the general interpretability of the decision-making process is particularly interesting for industrial applications because of reliability and acceptability issues.

2. Scalar feature extraction and risk of information loss. Classical machine learning algorithms such as DT, SVM and ANN take as input a vector of scalar features. This approach may be well suited to the analysis of steadystate disturbances. However, the time dependence of electrical waveforms is important for analyzing short-duration disturbances such as voltage sags.

The information related to the underlying cause is encoded through the entire duration of the event, and extracting scalar features involves a risk of information loss [START_REF] Susto | Chapter 9 -Time-Series Classification Methods: Review and Applications to Power Systems Data[END_REF]. There is also a higher risk of error in the calculation of scalar features if their extraction depends on a precise segmentation of transients and steady-states [START_REF] Caujolle | Impact of Waveform Segmentation Accuracy on Disturbance Recognition Reliability[END_REF].

For this reason, a classification approach based on time series seems more relevant for an efficient analysis of voltage sags. From the methods presented above, only [START_REF] Sha | Identification Method for Voltage Sags Based on K-means-Singular Value Decomposition and Least Squares Support Vector Machine[END_REF][START_REF] Yikun | Classification Method of Voltage Sag Sources Based on Sequential Trajectory Feature Learning Algorithm[END_REF] and the deep learning-based methods [START_REF] Li | Self-Supervised Voltage Sag Source Identification Method Based on CNN[END_REF][START_REF] Zheng | Recognition method of voltage sag causes based on two-dimensional transform and deep learning hybrid model[END_REF][START_REF] Zheng | Recognition method of voltage sag causes based on Bi-LSTM[END_REF][START_REF] Wang | Method of Voltage Sag Causes Based on Bidirectional LSTM and Attention Mechanism[END_REF] proposed a time series classification approach, without the need to calculate and extract scalar features.

Finally, it should be noted that all but one [START_REF] Liu | Identification and Location of Voltage Sag Sources Based on Multi-Label Random Forest[END_REF] of the previous literature references, target all the three main events responsible for voltage sags, regardless of their relative location (upstream or downstream). Nevertheless, industrial customers need to know whether the voltage drop is due to an event in the upstream distribution grid or the downstream industrial network in order to implement appropriate countermeasures. Indeed, if voltage sags are originated at the industrial site (downstream), corrective solutions can be more easily implemented to avoid future sags. However if the sags are originated at the level of distribution grid (upstream), corrective solutions are more difficult and very expensive to implement. In this case, industrial customers will prefer targeted mitigation solutions to protect sensitive equipment on their sites.

Conclusion

Power quality analysis is divided in three categories: detection, location and classification. Each category has a different goal, but they share common methods. We made a focus on the classification task and we presented a general five-stage scheme: data acquisition, pre-processing, feature extraction, feature selection and feature analysis. We also presented the most relevant techniques proposed in the literature for the classification of PQ disturbances.

The classification of voltage sag causes follows the same five-stage scheme and shares similar methods. Although the global classification accuracy of the proposals is globally high, we pointed out some important limitations that should be addressed in the development of our proposal:

• The large majority of the state of the art methodologies use statistical classifiers, which are highly data-driven. Thus, a large amount of data is required for the training of these type of classifiers. Access to these data is difficult, and therefore it should be limited.

• The generalization capabilities of these algorithms have not been evaluated, since all the methods have been trained partially or entirely with data obtained from the same source as the test data. This approach is not always feasible from an industrial implementation point of view, since real data is difficult to obtain and customers are not always capable of training the algorithms themselves. We should aim to develop an automated, scalable and generalizable methodology to be applied in different industrial sites. It should ideally be trained entirely on synthetic data prior to its deployment.

• Most of the extracted features are not easily interpretable from an electrical point of view. This can make the error analysis and trouble-shooting stage difficult to perform. Interpretable algorithms are preferable for acceptability and reliability issues. We should aim for an algorithm that integrates expert knowledge by selecting physical meaningful features.

• The extraction of scalar features from non-stationary signals involves a risk of information loss, since the time dependence of electrical waveforms is not taken into account. Alternative feature extraction methods should be studied in order to take into account this aspect.

• Few methods proposed a combined classification and location approach. Information on the relative location of disturbances is key to industrial customers, and therefore it should also be addressed.

The literature already proposes a large variety of solutions that report high accuracy for the classification of voltage sag causes. Our goal is to propose a solution that achieves such standards, but that at the same time addresses the aforementioned limitations. The algorithm to be developed has commercial purposes, hence its implementation must be realistic from an industrial point of view.

From the above analysis, we decided to propose a solution based on the extraction of meaningful and interpretable features by integrating expert knowledge in the process. These features are multivariate time series signatures, which compared to scalar values, have a lower risk of information loss. Then, we propose a time series classification approach to classify the signatures. Although this method is data-driven, we demonstrate that the amount of data required for the development of the solution is significantly lower compared to other approaches in the literature. In addition, we will ensure that the required data is accessible through simulation and that the developed algorithm is able to perform effectively in different industrial networks.

Chapter 3

Classification of Voltage Sag Causes

Introduction

In Chapter 2 we presented the most relevant methods and techniques applied to power quality analysis and the classification of voltage sag causes. Most of the proposals in the literature are based on the analysis of scalar features, with only a few algorithms making use of time series. Some of the major drawbacks of these methods include the risk of information loss linked to the extraction of scalar features, the low interpretability of the decision-making process and the high amount of data required for the training.

Therefore, inspired by pattern recognition, we propose a voltage sag classification algorithm based on the recognition of multivariate time series signatures. The signatures are obtained through a feature extraction stage mainly based on the Fortescue Transform (symmetrical components), which presents many advantages in terms of electrical interpretability. Although this transform has been applied in the literature for the characterization [START_REF] Bollen | Different methods for classification of three-phase unbalanced voltage dips due to faults[END_REF][START_REF] Oubrahim | Disturbances Classification Based on a Model Order Selection Method for Power Quality Monitoring[END_REF] and location [START_REF] Polajžer | Instantaneous positivesequence current applied for detecting voltage sag sources[END_REF][START_REF] Mohammadi | Employing instantaneous positive sequence symmetrical components for voltage sag source relative location[END_REF] of voltage sags, to the best of our knowledge, the obtained patterns have not been exploited for the classification of voltage sag causes. The signatures are classified using an approach based on the distance between an unknown signature and the labeled signatures in a previously built reference database. This approach has the advantage of requiring a reduced amount of data for its implementation. This chapter will describe the characteristics of our proposal. Each stage of the algorithm is presented, but the results and performance analysis will be further detailed in Chapter 4.

The first stage is the data acquisition of three-phase voltage and current waveforms. There are two sources of data available: synthetic data from numerical simulations and field measurements. In both cases, first, the data are pre-processed to ensure completeness, uniform length, and a correct sampling rate for all waveforms. Then, the Short-Time Fourier Transform (STFT), Fortescue transform, and Minmax normalization are applied to the waveforms in the feature extraction stage, providing a 4-dimension time series signature. At the end of this stage, the signatures of the labeled events are used for the constitution of the reference database.

The final stage is the feature analysis, which consists of calculating the distance between a new signature and the ones in the reference database. To calculate the distance between two signatures, we first align them in time and space using the Dynamic Time Warping (DTW) algorithm. The distance of the new signature to the different classes in the database is estimated, and the predicted label is given according to the closest class to the signature to classify. Finally, a confidence index associated with the predicted label is computed. The final output of the algorithm consists of the predicted voltage sag source and a confidence score.

Three main events are responsible for voltage sags: line faults, transformer energizing, and induction motor direct startup. These events can be generated in the distribution network (upstream) or the industrial network (downstream). While line faults and transformer energizing have been considered at both upstream and downstream locations, the induction motor direct startup has only been considered downstream due to its limited impact on the voltage of parallel feeders in LV sites. We also distinguish between balanced and unbalanced faults, depending on the fault causing the sag (see Chapter 1). As a result, we define seven classes for our classification problem, as detailed in Table 3.1.

Data acquisition

The final goal of this research work is to obtain an algorithm that will be applied to the data collected from different industrial sites. Each industrial site will be equipped with a single monitoring device on the LV side of the MV/LV main transformer. The algorithm could be deployed in two ways: in a centralized server where the data of several industrial sites would be analyzed, or directly onsite where the data are collected and processed locally (decentralized approach). Fig. 3.2 illustrates both approaches. The algorithm should ideally be developed using only real field data. However, as discussed in Chapter 2, a significant drawback of this type of data is the difficulty of collecting a sufficient number of PQ events because they are power system anomalies and, therefore, rare. Another limiting aspect of field data is the labeling of events, which requires expert analysis. In addition, the human factor must be considered, as there is a risk of error in the labeling process.

Therefore, it is preferable to rely on synthetic data to create a large and representative dataset. It allows us to obtain large amounts of labeled and controlled measurements, even for events that are so rare that they have not been registered yet or for which we have few records. Moreover, since we can control different parameters, we can evaluate the algorithm's robustness regarding noise levels or frequency variations.

Real field data

Field data were acquired at three industrial sites over a period of 1.5 years, each from a different sector: metal equipment manufacturing, agri-food industry, and chemical manufacturing. These LV sites are fed from the 21 kV distribution network via a 21 kV/400 V Dyn transformer. It is important to note that no prior knowledge of the industrial network topology or equipment was used to develop this algorithm, as this information is unavailable.

The voltage sag detection algorithm, which is based on threshold analysis, is built into the monitoring device at factory. Once a voltage sag is detected, the voltage and current waveforms for 16 cycles are recorded with a sampling frequency of 12.8 kHz1 . As presented in Table 3.2, the field dataset consists of 385 voltage sags, belonging to 5 available classes out of the 7 considered: Upstream balanced faults (A1), upstream unbalanced faults (A2), downstream unbalanced faults (B2), upstream transformer energizing (C1) and downstream motor starting (D). The classes corresponding to downstream balanced faults (B1) and downstream transformer energization (C2) were not recorded during this monitoring period. The data labeling has been validated by two power quality experts. Note that the spread of the events across classes and sites is not uniform. Indeed, upstream unbalanced faults are the most common source of voltage sags, followed by upstream balanced faults (three-phase faults). However, their occurrence depends on the characteristics of the distribution network feeding the site, and the three monitored sites are geographically located in different regions of the country (different feeder lines). Similarly, the frequency of sags caused by upstream transformer energizing will depend on the distribution grid characteristics.

On the other hand, downstream faults are less frequent. Voltage sags caused by the direct startup of downstream induction motors are due to oversized motors without starting devices and depend on how the site has been sized. Finally, isolation transformers are not always present in LV sites and should be oversized to cause voltage sags.

Synthetic data generation

Characteristics of the simulated industrial grid

The simulated case study is designed using the EMTP-RV (ElectroMagnetic Transients Program) software. The model comprises two sub-networks: the industrial network and the distribution network, as displayed in Fig. 3.3. It should be noted that the monitoring device is located downstream of the main MV/LV transformer. Thus, events occurring in the distribution network are qualified as upstream (regarding the monitoring point), and those occurring in the industrial network as downstream.

Industrial network

For our study, the nominal voltage of the industrial network is 400 V, and the frequency is 50 Hz. The site is connected to a Dy11n 21 kV/400 V transformer of 400kVA. We have included the loads described in Table 3.3 to model the behavior of a small but representative LV industrial site. The loads are not all simultaneously connected. For instance, the large motors and transformers are connected only for the motor startup and the transformer energizing scenarios.

The total power consumption of the modeled industrial site is lower than 380 kVA. The main 400 kVA transformer is not optimally sized since we want to emulate a critical state where the power grid is susceptible to experiencing voltage sags while maintaining a total load just under the nominal transformer capacity. Finally, we note that some of the devices are sources of permanent disturbances, frequently present in industrial networks:

• Harmonics: Variable speed drivers and 6-pulse rectifiers generate current harmonics, mostly of order 5 and 7.

• Unbalance: The presence of single-phase loads simulates the injection of a moderate level of unbalance.

Distribution network

The distribution or upstream network has a voltage level of 21 kV. The main busbar's nominal power varies from 10 to 100 MVA, and the ratio of the reactance to the system's resistance is 0.1. The length of the MV lines varies between 0 and 30 km. Two feeder lines 2 connect the 21 kV busbar to the main industrial network and to a secondary site. The latter is connected to the grid through two Dyn11 transformers with power ratings ranging from 500 to 1250 kVA. The transformers are connected in parallel, each feeding a 300 kW load.

Generation of voltage sags

The generated voltage sags vary between 10% and 98% (residual voltage). Although a voltage sag is defined as the reduction of RMS voltage under 90% of its rated value, shallow voltage drops are also included in the dataset. Their identification can be more complex because the voltage and current variations can be very small. A description of the parameters used for the generation of voltage sags is given below:

• Line faults Different line fault types are generated: single line-to-ground faults (L-G), line-to-line faults (LL), double line-to-ground faults (LL-G), and three-phase faults (LLL) and three-phase-to-ground faults (LLL-G). The magnitude of the voltage sag mainly depends on the fault type, the distance to the fault (line's length), and the value of the ground fault resistance. The duration of the sag, in this case, depends on the action of protective devices. The fault clearance time varies between 50, 100, and 200 ms to obtain sags of varying duration. We also varied the distance between the monitoring and fault injection points to obtain a significant variation in voltage dip amplitudes. For upstream faults, we varied the distance between 0 and 30 km, and the value of the ground fault resistance between 0 and 40 Ω. For downstream faults, we varied the distance between 0 and 20 m and the fault resistance between 0.1 and 10 Ω. This type of event is generated at three locations: the MV feeder of the secondary site (upstream), the 400 V busbar of the main industrial network (downstream), and the secondary side of an isolation transformer (downstream).

• Transformer energizing

The magnitude of the voltage sag caused by this event depends on the transformer's power and its initial core's flux (different from zero when the core has not been entirely demagnetized before re-powering). The duration of the sag depends on the transformer's characteristics. These events were generated using transformers with power ratings of 500 and 1250 kVA for upstream transformers and 100 to 250 kVA for downstream transformers. We varied the initial core's flux and the energizing starting time within a cycle duration (0 ms, 5ms, 10ms, and 18ms). This type of event is generated at two locations: the secondary side (upstream) and the 400 V busbar of the industrial network (downstream).

• Induction motor direct startup

The magnitude of the voltage sag caused by a motor startup depends on its power, torque, and total inertia. These parameters also determine the duration of the event and the sag amplitude. We simulated various profiles of motor direct startups using induction motors with power ratings varying from 22 to 110 kW. The load torque and inertia are also varied according to the motor's size. This type of event is generated by directly connecting the induction motor to the 400 V busbar of the industrial network (downstream). It should be noted that the presence of variable speed drivers, three-phase rectifiers, and unbalanced loads simulate harmonic distortion and unbalance.

We have generated a synthetic dataset of 100 events for each of the seven defined classes (i.e., 700 events). The sampling frequency of the synthetic data is set at 12.8kHz to match the settings of the actual monitoring device. However, we will see later that this value can be reduced up to 400 Hz without degrading the algorithm's performance.

Pre-processing

At this stage, we perform simple operations on the voltage and current waveforms to ensure the uniformity of the data. First, the sampling frequency of all the voltage and current waveforms is recalculated, and if needed, they are resampled to match the reference value (12.8 kHz). The waveforms are also trimmed to obtain a uniform length.

Then, the sag's completeness is verified since the monitoring device may sometimes record incomplete voltage sags. Ideally, there is at least one healthy cycle previous to the sag. However, the proposed solution can tolerate up to a single healthy half-cycle before the sag. Finally, to avoid the "loss" of a cycle when calculating the Short-Time Fourier transform in the feature extraction stage, we artificially add one pre-sag period by copying the first cycle in voltage and current. The reason for this "lost cycle" is due to the moving window in the STFT since it uses the first cycle to calculate the first data point of the transformed signal.

Feature extraction: Multivariate Time Series Signature

The main objective of this step is to bring out features from the voltage and current waveforms that are common to events from the same class but different from those in other classes. Transforming the raw waveforms into simpler but meaningful components helps and simplifies the classification task. Indeed, using raw waveforms or even RMS curves directly would require classification algorithms of a higher level of complexity, requiring larger amounts of training data. In addition, the decision-making process would be more challenging to interpret.

In this stage, the pre-processed voltage and current waveforms are transformed into a 4-dimension time series signature. This transformation combines the Short-Time Fourier Transform (STFT) and the Fortescue transform. Finally, a min-max normalization is applied to the extracted signatures since the classification is based on the shape of the signatures. Each step is detailed in this section.

Short-Time Fourier Transform

STFT is first used to decompose the signal into its harmonic components, extracting only the harmonics of interest and avoiding those affected by the industrial loads. For instance, harmonics of orders 5 and 7 are affected by 6-pulse rectifiers, present in variable speed drivers and three-phase rectifiers. Their elimination increases the generalization capability of the algorithm despite the different load types present in the industrial network. This also allows us to reduce the number of simulated scenarios for the "training" stage compared to other methods in the literature that take as input raw waveforms, including all the harmonic content.

The fundamental frequency (first harmonic, f = 50Hz) contains prime information for the analysis of voltage sags because sags are low-frequency disturbances. The presence of even harmonics due to transformer energizing is also helpful information for its identification. Thus, the = second harmonic (f = 100Hz) is also retained. Finally, the first two harmonics will be used for the analysis of voltage sags.

The second reason for implementing STFT is that the Fortescue transform (which is applied next) is defined in the complex domain. Thus, at the end of this step, the voltage and current waveforms are converted into two complex harmonic phasors.

The Short-Time Fourier Transform (STFT) [START_REF] Kehtarnavaz | CHAPTER 7 -Frequency Domain Processing[END_REF] is a technique used to analyze the frequency content of a non-stationary signal. It is performed by applying the Discrete Fourier transform (DFT) to the discrete time signal through a sliding window of length W L . The window overlap between the signal segments compensates for the signal attenuation at the window edges. The DFT of each segment is stored in a matrix containing the magnitude and angle for each data point in time and frequency.

The matrix is defined in equations (3.1) and (3.2).

ST F T

(x) = [X 0 (f ), X 1 (f ), ..., X T (f )] (3.1) X m (f ) = ∞ n=-∞ x(n)g(n -sW H )e -j2πf n (3.2)
Where m is the window number, n is the length of the discrete signal x(n), X(f ) is the DFT of the windowed signal centered at time sH, W H is the hop length between successive windows3 and g(n) is the window function.

Fortescue transform

The choice of the Fortescue transform compared to other transformations such as the Wavelet transform (WT), S-transform (ST), Hilbert-Huang transform (HHT), and even the Clarke transform is due to its electrical interpretability. Indeed, electrical equipment can be modeled using equivalent circuits based on symmetrical components [START_REF] Das | Understanding Symmetrical Components for Power System Modeling[END_REF]. For instance, Milanovic et al. [START_REF] Milanovic | The Influence of Induction Motors on Voltage Sag Propagation-Part I: Accounting for the Change in Sag Characteristics[END_REF] presented an analytical approach for analyzing the interaction between induction motors and voltage sags. The analytical equations derived from the equivalent circuits based on symmetrical components proved to accurately reproduce the behavior of the corresponding transient simulation model. Moreover, the voltage sag's 4-D signature that we propose for class A2 (upstream unbalanced faults) is found and analyzed in detail by the authors in this reference.

Fortescue transform [START_REF] Fortescue | Method of Symmetrical Co-Ordinates Applied to the Solution of Polyphase Networks[END_REF] is a linear transformation used to analyze unbalanced three-phase power systems. It transforms an unbalanced set of three phasors (X A , X B , X C ) into a balanced set of three symmetrical components: the positive (X 1 ), negative (X 2 ) and zero-sequence (X 0 ), as illustrated in Fig. 3

  X A X B X C   =   X a0 X b0 X c0   +   X a1 X b1 X c1   +   X a2 X b2 X c2   (3.3)
The three symmetrical systems are perfectly balanced, with an angle equal to 2 3 π between each phase. Therefore, each phase in the system can be written in function of the phasor rotation operator α = e j 2 3 π and X a0 , X a1 , X a2 , which will be simply referred as X 0 , X 1 , X 2 . Note that the positive sequence can also be noted as X + and the negative sequence as X -. The Fortescue transform of a three-phase system for a given harmonic h is defined in (3.4).

  X 0 X 1 X 2   {h} = 1 3   1 1 1 1 α α 2 1 α 2 α     X A X B X C   {h} (3.4)
The Instantaneous Symmetrical Components or ISCs (X 0{h1} (t), X +{h2} (t), etc.) of the harmonics extracted from the voltage and current waveforms are calculated using the instantaneous complex values previously determined by STFT, as illustrated in Fig. 3.9. The positive-sequence component of the fundamental frequency represents the actual voltage and current being provided to the load. In a perfectly balanced system, the negative and zero-sequence components are equal to zero. The zerosequence is directly related to the grounding system and transformer winding connections. Thus, an increase of this component is visible for line-to-ground faults (L-G, LL-G, LLL-G) but not for line-to-line faults (LL or LLL). This is true for downstream faults, however, for upstream disturbances, the zero-sequence component is filtered [START_REF] Aung | The Influence of Transformer Winding Connections on the Propagation of Voltage Sags[END_REF] since MV/LV transformers of industrial sites usually have a Dy winding connection. No zero-sequence current is generated because the induction motors' windings are connected in delta (D) or ungrounded wye (Y ). Thus, the zero-sequence voltage is not influenced by the induction motor either [START_REF] Yaleinkaya | Characterization of voltage sags in industrial distribution systems[END_REF].

For these reasons, four ISCs are selected to compose a 4-dimensional time series signature:

(a) Voltage positive-sequence, harmonic 1 (b) Voltage negative-sequence, harmonic 1 (c) Voltage positive-sequence, harmonic 2 (d) Current positive-sequence, harmonic 1

The characteristics of each ISC depending on are described in more detail in Section 3.5.4.

Min-max normalization

The selected ISCs constitute a 4-dimension time series signature. We use the shape of each ISC as the main discriminant characteristic between classes. Thus, the amplitude difference among events from the same class must be reduced since it is not a discriminant criterion. This should be achieved without degrading the original shape. Therefore, the most suitable normalization method is minmax, which rescales the entire time series proportionally within a defined value range. It also preserves the shape characteristics of the signature contrary to zscore normalization, whose goal is to obtain a zero mean and a unitary standard deviation.

For this reason, we apply a min-max normalization to each ISC to perform a shape-based time series classification. Minmax normalization with a re-scaling between [l a , l b ] is defined in (3.5) with X being the ISC to be normalized, l a = -0.5 and l b = 0.5.

X [la,l b ] = X -min(X) max(X) -min(X) * (l b -l a ) + l a (3.5)
The ISCs are first re-scaled, then zero-centered by subtracting the first point value from the rest of the sequence. Each time series is defined between [-1, 1] but with a spread equal to one at the end of this operation. Fig. 3.12 illustrates the feature extraction process and the obtained components at each stage, down to the four-dimensional signature considered.

To reduce the overall computation time, we apply a resampling operation of 1/32 to each ISC. This is possible because, as mentioned previously, the minimum sampling rate for the waveforms can be reduced up to 400 Hz without degrading the algorithm's performance because the Shannon theorem F s ≤ 2 • f max is still respected, with f max = 100Hz (second harmonic). As an example, Fig. 3.10 illustrates the third ISC of the signatures obtained from downsampled waveforms. The error due to the downsampling can be estimated by calculating the Normalized Root Mean Squared Error (NRMSE) between the obtained signatures. For this, we used the 700 events in the synthetic dataset. The NRMSE in two dimensions (L, D) is defined in (3.6). Where L is the length of each ISC, D the number of dimensions, s the original signature at 12.8 kHz and s ′ the signature obtained from the downsampled waveforms. Because we applied a minmax normalization, the spread of each ISC is |l b -l a | = 1. Fig. 3.11 summarizes the RMSE between the reference signatures obtained from the original waveforms at 12.8 kHz, and those obtained from the downsampled waveforms at 6.4 kHz, 400 Hz, and 200 Hz. The results show that a sampling rate up to 400 Hz is possible with an error of 4.02% in the feature extraction process. The four-dimension signatures of voltage sags from field data records are illustrated: Fig. 3.13 corresponds to an upstream unbalanced fault, Fig. 3.14 corresponds to an upstream transformer energizing, and Fig. 3.15 corresponds to a downstream motor direct startup. It should be noted that the four components of the signature or ISC are correlated in time. The differences between the three events are more visible with the signatures than with the raw waveforms. The simplicity of these signatures is an advantage to the classification stage. Moreover, Fig. 3.16 illustrates the advantages of the feature extraction process for bringing out relevant and common characteristics between events of the same class (upstream unbalanced fault). Although there are notable similarities in the waveforms and RMS curves, the differences in amplitude, the affected phases, and the harmonic content are not negligible. Using this type of data as input to the classifier would require the implementation of complex algorithms. However, these differences are minimized when their signatures are compared, with only a significant difference in the fault duration. The signatures are particularly invariant to the variation of the sag's amplitude. These characteristics will allow the use of a simple but effective approach for the classification step.

Electrical interpretation of the voltage sag signatures

The four-dimensional signatures obtained at the end of the pre-processing stage are interpretable from an electrical point of view. We detail in the following section the main characteristics of each ISC. The full list of signatures corresponding to each class is illustrated in Appendix A.

1. Voltage positive-sequence harmonic 1 (Fig. 3.17). This component represents the voltage evolution during the sag, with very similar characteristics to the RMS three-phase voltage curves. The rapid drop is directly related to the occurrence of the voltage sag. For example, the quasi-square shape is characteristic of fault-induced voltage sags. Recovery is rapid after the fault is cleared, as shown in Fig. 3.17a. Sags caused by starting motors and energized transformers have, in contrast, a progressive recovery with characteristics similar to their RMS voltage curves, as illustrated in Fig. 3.17b and 3.17c. This first component helps distinguish between drops due to a fault, the energizing of a transformer, or the direct starting of a motor. However, this component alone does not provide sufficient information to determine whether the fault is balanced or unbalanced, or where it is located. We use this component to determine the relative location of the event (upstream or downstream of the monitoring device) [START_REF] Tayjasanant | A Resistance Sign-Based Method for Voltage Sag Source Detection[END_REF][START_REF] Polajžer | Instantaneous positivesequence current applied for detecting voltage sag sources[END_REF][START_REF] Mohammadi | Employing instantaneous positive sequence symmetrical components for voltage sag source relative location[END_REF]. This component increases and remains constant for the duration of the voltage sag for downstream events. It then decreases proportionally to the voltage sag recovery curve. Fig. 3.20a, 3.20b and 3.20c illustrate this for the cases of a downstream fault, a downstream transformer energization and a downstream direct motor start. For upstream events, on the contrary, it rapidly decreases when the event starts. This phenomenon can be explained by the energy sink analogy presented in [START_REF] Parsons | A direction finder for power quality disturbances based upon disturbance power and energy[END_REF], where events such as faults or load connections consume high amounts of current and energy. Then, the current recovers and stabilizes, particularly in the case of fault-induced sag. When the voltage is restored at the end of the sag, the current will reach a peak value higher than its nominal value. The magnitude of the peak depends on the duration and severity of the voltage drop and the connected loads (motor power and inertia, DC bus capacity, etc.) Fig. 3.20d and 3.20e illustrate the current component corresponding to an upstream fault and an upstream transformer energizing respectively. As previously mentioned, several methods in the literature use classical machine learning algorithms that mainly exploit scalar features. This approach is more suitable for the analysis of steady-state disturbances. However, the time dependence of electrical waveforms is important for analyzing short-duration disturbances such as voltage sags. Therefore, we propose a time series classification approach, exploiting the multivariate time series signatures obtained after the feature extraction process. Time series classification methods are usually based on: distance measures, shaplets, dictionaries, interval summarising or deep learning (neural networks). Among these methods, 1-Nearest Neighbor with Dynamic Time Warping (1NN-DTW) is a recommended benchmark due to its simplicity and hard-to-beat accuracy [START_REF] Bagnall | The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF][START_REF] Ruiz | The great multivariate time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF]. It also requires significantly less data than deep learning algorithms.

For these reasons, we have chosen the 1NN-DTW algorithm as our first approach for our task. A new voltage sag is then classified by calculating the distance between its signature and the ones in the reference database (distance-based classification). We first define a distance measure to compare multivariate time series signatures based on Dynamic Time Warping. However, one of the main drawbacks of 1NN-DTW is its low robustness to outliers since the prediction is defined based on a single neighbor. Therefore, we proposed an improved version of 1NN-DTW, which we refer to as Nearest Neighborhood classifier. It is based on the distance to an entire "neighborhood" or class, instead of a single neighbor. In addition, we present two methods for the distance-to-class calculation: a mean estimation method with bootstrapping and a centroid estimation method with soft-Dynamic Time Warping (soft-DTW). Finally, we define two confidence scores associated with the label prediction.

Distance between signatures

Voltage sags caused by the same event (i.e., belonging to the same class) have very similar signatures, as previously illustrated in Fig. 3.16. The signatures are invariant regarding sag's amplitude variations, thanks to Min-max normalization. Most notable differences are due to varying starting times and duration of the events, which result in time misalignments. Another possible source of misalignments between signatures is the imprecision in the detection algorithm within the monitoring device. Indeed, an imprecise detection can sometimes cause incomplete recordings, resulting in space misalignments. Therefore, we first perform a spatio-temporal alignment to compare the signatures' similarities properly.

Space alignment

The characteristics and parameters of the monitoring device may not always be perfectly calibrated. This can cause incomplete event recordings, affecting the signature's matching process. As presented in the pre-processing stage, a complete recorded voltage sag would ideally have one or more healthy cycles prior to the sag. However, thanks to the feature extraction process that we have proposed, it is still possible to evaluate the event if there is at least a 1/2 pre-sag cycle. In this case, an offset will be applied to correct the space misalignment.

We first check whether there is a risk of having less than one reference period due to an incomplete recorded event by calculating the difference between the average values of the first period and the last period for each ISC in the signature. If the difference δ is higher than a minimum threshold of 0.05 pu, an offset must be added to improve the alignment. A vector containing the possible offsets is defined as X of f = -sign(δ) * [0, 0.05, 0.1, ...δ] and the Euclidean distance is calculated between the ISC query plus the offset and the ISC reference. The minimal Euclidean distance gives the optimal offset, which is applied at the end of this stage.

Time alignment using Dynamic Time Warping Dynamic Time Warping for multivariate time series

Classical distance measures such as Euclidean distance can result in poor results because they do not consider the differences in speed and time delay between time series [START_REF] Bagnall | The great time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF][START_REF] Ruiz | The great multivariate time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF]. Therefore, we use the Dynamic Time Warping (DTW) algorithm proposed in [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] for handling time alignment differences between two univariate time series 4 . Fig. 3.21 illustrates the difference between the standard Euclidean distance and Dynamic Time Warping when matching two time series.

The alignment of two time series, r (reference) and q (query) of length L r and L q respectively, is given by the warping path ϕ = (ϕ(1), ..., ϕ(Lw)) of length L w . Where ϕ(l) = (ϕ r (l), ϕ q (l)) ∈ [1 : L r ] × [1 : L q ], for l ∈ [1 : L w ], satisfying the following conditions:

• Boundary condition: ϕ(1) = (1, 1) and ϕ(L w ) = (L r , L q )

• Monotonicity condition: ϕ r (l -1) ≤ ϕ r (l) and ϕ q (l -1) ≤ ϕ q (l)

• Step size condition: ϕ r (l) -ϕ r (l -1) ≤ 1 and ϕ q (l) -ϕ q (l -1) ≤ 1 Although the algorithm can handle time series of different lengths, the lengths of the compared signatures are uniform and will be referred to in the following as L. To calculate the optimal warping path that best aligns both times series, we first calculate a local distance matrix of size L × L, L being the length of the signature. Each element of the distance matrix is defined through a chosen distance measure, the Euclidean distance being usually the privileged choice. The equation calculating the elements a(i, j) of the distance matrix A for two univariate time series q (query) and r (reference) with indexes i, j ∈ [0, L] is given by equation (3.7). a(i, j) = (q(i) -r(j)) 2 (3.7)

The optimal warping path ϕ is found through the minimization of the cumulative cost E defined in (3.8) obtained from the distance matrix. It is given by two vectors of integers ϕ q (l), ϕ r (l) of same length L w (with L <= L w <= 2L), mapping the time axis of the query q to the reference r, and w(l) is an optional per-step weighting coefficient. These vectors indicate the time alignment to be applied to all the dimensions of both signatures q and r. Fig. 3.22 illustrates the principle of time alignment of two univariate time series. Extending this algorithm to multivariate time series can be achieved in two ways. We can independently calculate the DTW distance for each dimension and add the calculated distances (independent dynamic time warping DTW I ). We can also calculate the DTW distance across all the dimensions simultaneously (dependent dynamic time warping or DTW D ) [START_REF] Ruiz | The great multivariate time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF]. Both approaches are illustrated in Fig. 3.23. In our case, the four dimensions of each signature are time-correlated, and a single optimal time alignment warping path is calculated for all four dimensions. The equation to calculate the elements a(i, j) of the distance matrix A for two multivariate time series signatures q and r of dimension D = 4 with indexes i, j ∈ [0, L] becomes:

E(q, r) = min

a(i, j) = D=4 d=1 (q(i, d) -r(j, d)) 2 (3.9)
Finally, we define a normalized distance D(q, r) between the aligned query q and reference r, as defined in equation (3.10). The distance is normalized by the length L w . This normalization helps to obtain a similarity measure independent of the signature deformation since L w tends to increase when dilations or compressions are applied to the signature during time alignment.

D(q, r) = 1 L w D=4 d=1 Lw l=1 (q(ϕ q (l), d) -r(ϕ r (l), d)) 2 (3.10)

Local constraint optimization

DTW is a powerful tool for time series distance measurement. Nonetheless, it can lead to pathological alignments caused by the so-called "singularities", as illustrated in Fig. 3.24. A singularity is a data point from one time series that matches a large section of the second time series, leading to highly distorted sequences. This is an important issue for our problem since the classification is shape-based. Thus, the global shape of the signatures should be kept as close to its initial form as possible. Local or global constraints can be applied to limit this phenomenon. These constraints are applied by modifying the per-step weights w(l) used to calculate the cumulative cost.

Global constraints such as the Sakoe-Chiba band [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] or the Itakura parallelogram [START_REF] Itakura | Minimum prediction residual principle applied to speech recognition[END_REF] limit the distance of the warping function to the main diagonal, as shown in Fig. 3.25. However, these constraints are not suitable for our problem as they may prevent an optimal alignment of two similar events but with very different durations, thus limiting the time dilation (stretching) or compression capability of the DTW algorithm. Local constraints include step patterns, which are more flexible but still can limit severe signal distortion.

Step patterns describe the per-step weighting coefficients w(l) to include in the cumulative cost equation in (3.8). Fig. 3.26 illustrates three well-known step patterns: symmetric2, symmetric1 and asymmetric [START_REF] Giorgino | Computing and Visualizing Dynamic Time Warping Alignments in R : The dtw Package[END_REF]. The numbers on the transitions indicate the multiplicative per-step weight. Note that the asymmetric step pattern is the only one that does not respect the step size condition, but it is still considered a valid step pattern in the literature. Sym-metric2 is the default choice in various DTW implementations [START_REF] Giorgino | Computing and Visualizing Dynamic Time Warping Alignments in R : The dtw Package[END_REF]. The selection of the best step pattern requires the definition of a coefficient that allows us to compare the time alignment performance of the step patterns. However, classical error metrics such as mean absolute percentage error (MAPE) can not be directly applied for this purpose since it requires knowing the optimal warping path or ground truth for a given time alignment between two signatures, which is impossible to know in practice.

Therefore, for our application, we look for the best trade-off between distance and signature distortion. The best step pattern should allow minimal distance between signatures of the same class and maximal distance between signatures of different classes, while assuring a minimal distortion of the time series in both cases.

The distance D(q, r) between two signatures is given by equation (3.10). We defined a distortion rate Z (%) in equation (3.11) based on a similar coefficient presented in [START_REF] Keogh | Derivative Dynamic Time Warping[END_REF], with L w being the warping path length, and L the length of the original signature. Since the two compared signatures have the same length, a perfect time-alignment would correspond to a diagonal warping path of length L w = L and Z (%) = 0. The maximum warping path length being L w = 2L, the maximum value of Z (%) = 100.

Z (%) = L w -L L (3.11)
We have compared these three step patterns using the distortion rate Z, the distance D between signatures belonging to the same class (intra-class distance), and the distance between signatures of different classes (inter-class distance). From the results summarized in Table 3.4 5 , we can conclude that symmetric1 exhibits the best trade-off between minimal intra-class distance, maximal inter-class distance, and minimal distortion rate. The results are better illustrated in Figs. 3.27 to 3.29. For instance, Fig. 3.27 illustrates the time alignment result using the sym-metric2 and symmetric1 step patterns between two signatures belonging to the same class (intra-class). We note that the time alignment result is correct in both cases. This is reflected in Table 3.4, where the intra-class distances of symmetric2 and symmetric1 are very close (1.41 and 1.85 respectively). However, the warping path is smoother and closer to a diagonal for symmetric1 compared to symmet-ric2, resulting in a much lower distortion rate (19.79% for symmetric1, compared to 71.83% for symmetric2 ). The time alignment in Fig. 3.28 is performed between signatures belonging to different classes (inter-class). When using the symmetric1 step pattern, the DT W D algorithm attempts to align both signatures while preserving the original shape, resulting in a low distortion rate. This is not the case for symmetric2, where the query signature is significantly distorted to fit the reference signature. This high distortion is due to the presence of multiple singularity points, corresponding to a warping path with long vertical and horizontal segments. Moreover, such distortion within the signatures can lead to false low distances D between signatures of different classes, resulting in misclassification errors due to poor time alignment. The results in Table 3.4 illustrate this issue: the symmetric2 step pattern has a higher distortion rate of 90.11% and a lower inter-class distance of 6.26, compared to symmetric1 with a distortion rate of 20.59% and an inter-class distance of 9.48. Finally, Fig. 3.29 shows the results when comparing the asymmetric and sym-metric1 step patterns. The asymmetric step pattern fails to align both signatures since it has a reduced dilation capacity. On the other hand, the symmetric1 step pattern allows a more flexible stretching of the query signature resulting in better results. This is also reflected in Table 3.4, with an inter-class distance of 2.26 for asymmetric and 1.85 for symmetric1.

In summary, the symmetric1 step pattern favors oblique steps over horizontal or vertical ones. This characteristic limits distortion compared to the symmet-ric2 step pattern, which considers a vertical plus horizontal step equivalent to an oblique step. The symmetric1 step pattern also allows a higher degree of dilation or compression compared to the asymmetric step pattern since the latter imposes a single match point for each point in a time series. Although DTW is a well-known algorithm in fields such as speech processing and time series analysis, it has only been implemented, to our knowledge, once for the classification of types of PQ disturbances by Youssef et al. in [START_REF] Youssef | Disturbance Classification Utilizing Dynamic Time Warping Classifier[END_REF]. The authors used Fast Fourier transform (FFT), and Walsh transform for feature extraction and a combination of Vector Quantization (VQ), fast matching technique, and DTW for the identification of six classes of numerically simulated disturbances, including the following classes: no disturbance, voltage swell, voltage sag, harmonics, voltage flicker, and oscillatory transient. The success rate was 97%, and the fast matching technique helped reduce the overall computation time by 66.6%. However, a significant drawback is the number of disturbances needed to constitute the training dataset, with 200 events per class.

Distance-based classification

We present in this section two distance-based classification methods: 1-Nearest Neighbor (1NN) classifier and Nearest Neighborhood classifier. Both classifiers use the distance D(q, r) defined in equation (3.10), based on the DTW. The 1NN-DTW classifier is first explored as a benchmark. Then, to overcome the limitations of 1NN-DTW, we propose an improved version of this classifier based on the distance to an entire "neighborhood" or class instead of a single neighbor. Two methods for estimating the distance to a given class are also presented.

1-Nearest Neighbor classifier

The 1-Nearest Neighbor with Dynamic Time Warping (1NN-DTW) classifier was first selected since it is a recommended benchmark for time series classification due to its simplicity and hard-to-beat accuracy [START_REF] Ruiz | The great multivariate time series classification bake off: a review and experimental evaluation of recent algorithmic advances[END_REF]. However, a requirement for better results when using distance-based classifiers is to have a good class separability. This aspect will be studied in further detail in Chapter 4.

1NN-DTW is based on the k-Nearest Neighbors algorithm, κ representing the number of closest neighbors (for 1NN, κ = 1). It works by calculating the distances between a query (new event to classify) and the references in the training dataset. It selects the closest reference to the query, and the selected reference's label gives the classification result. Considering that it is a distance-based classifier, the choice of the distance measure is essential. To properly calculate the distance between two multivariate time series signatures, we use the distance D(q, r) defined in (3.10), based on the dependent Dynamic Time Warping algorithm (DT W D ). Fig. 3.31 illustrates the principle of classification using 1-Nearest Neighbors.

The efficiency and simplicity of this solution makes it a good first approach. However, there are some limitations because the classification is only based on a single neighbor. Indeed, the robustness of the algorithm is not guaranteed, since the presence of one or more outliers in the training or reference database may cause misclassification errors. In addition, k-NN does not provide a probability or confidence index associated with the prediction. The closest method to calculate a similar score is a majority voting system. For this, the number of closest neighbors κ must be set to a value higher than 1. The result would be the number of neighbors of each class, divided by the total of κ neighbors. For instance, take κ = 5, with three neighbors belonging to class 1 and the other two neighbors belonging to class 2 and class 3. The "probability" of the new event belonging to class 1 would be 3/5 (or 60%), 1/5 (or 20%) for class 2 and 1/5 (or 20%) for class 3. This estimate of the probability of membership is neither accurate nor sufficiently robust enough.

Finally, there is also a limitation in terms of computation speed. Searching for the closest neighbor requires calculating the distance to all the events in the reference database. Real-time calculation is not a primary constraint, given that the mitigation solutions can only be deployed after analyzing a certain amount of historical data. However, a reduced computation time is still desirable. This can be particularly useful if the classification algorithm is to be implemented in a decentralized platform with limited hardware resources.

For these reasons, we propose an improved solution based on this first approach to overcome these limitations.

Nearest Neighborhood classifier

The sensitivity of the 1NN-DTW classifier to outliers can be reduced if we estimate the distance to an entire class or "neighborhood". Therefore, to increase the robustness of the 1NN approach, we consider the totality of the signature references in the database. Thus, a new signature is classified by calculating the distance to entire groups of signatures belonging to the same class. The new signature is then labeled according to the closest class. We propose two methods for calculating the distance to a particular class of a new sag in the following sections: a mean distance estimation with bootstrapping and a distance calculation with centroid estimation.

Mean distance estimation with bootstrapping

The most intuitive method for calculating the distance of a new voltage sag signature q * to a class C k of size N k , is to calculate the mean distance to all the reference signatures r k in the class. However, the mean is not a robust estimator for outliers. Thus, to obtain a more robust estimator of the mean distance to each class d k (q * ), we define dk (q * ) using a bootstrapping approach as described in (3.12), where B is the total number of sub-samples X D,b extracted from the population of distances X D = {D(q * , r k,n )} with n = 1, 2, ..N k , and dk,b is the mean of X D,b as defined in (3.13). It can be noted that this is a resampling with replacement technique. Therefore the size of the sub-samples B can be equal to the number of signatures per class N k , which is the case here (B = N k ). Fig. 3.32 illustrates the distance calculation strategy using bootstrapping.

d k (q * ) ≈ dk (q * ) = 1 B B b=1 dk,b (3.12) dk,b (q * ) = 1 N N n=1 D(q * , r k,n ) (3.13)
Once the distances d k (q * ) between the new voltage sag signature q * to each class C k are calculated, we obtain a vector containing the distances to the K classes. The closest class gives the label y * assigned to the event, such as: With this approach, the sensitivity to outliers can be effectively reduced. However, the total computational cost remains very high. The cost due to the distance calculation to each element in the database is O(N k • K • D • L 2 ), N k being the number of signatures per class, K the number of classes, D the dimensions of the signature and L the length of the time series. Therefore, a second method for distance-to-class calculation is explored.

y * = arg min k∈{1,...,K} ([d 1 (q * ), d 2 (q * ), d k (q * )..., d K (q * )]) (3.14)

Centroid estimation with soft-DTW

The main disadvantage of the mean distance estimation with bootstrapping is the computational cost, since the method requires the calculation of the distances to all the signatures in the database. The distance of a new voltage sag signature q * to each class can also be obtained by reducing the class to a single representative signature. This signature can be assimilated as the centroid of the class. The classification of a new voltage sag would be performed by calculating the distance of the new signature to each centroid in the database. Then, the new signature is labeled according to the closest centroid.

In terms of computational cost when evaluating a new voltage sag, the cost of 1NN-DTW and Nearest Neighborhood using mean distance with bootstrapping is

O(N k • K • D • L 2 )
, N k being the number of signatures per class, K the number of classes, D the dimensions of the time series and L the time series length. The new approach reduces the computational cost by N k because each class is now represented by a single centroid. The distance d k (q * ) between a new signature q * and a class C k is given by the distance to its centroid.

Nevertheless, determining the centroid of a multivariate time series cluster is not trivial. Calculating the barycenter of a set of time series requires time alignment, as achieved with Dynamic Time Warping (DTW). For this reason, we have implemented a state-of-the-art variant of DTW for barycenter calculation: soft-Dynamic Time Warping (soft-DTW). Fig. 3.33 illustrates the Nearest Neighborhood with a centroid estimation approach using soft-DTW. It is to be noted that the centroid estimation is performed once and in "off-line" mode, prior to the implementation of the algorithm during "on-line" operation. Soft-DTW [START_REF] Cuturi | Soft-DTW: a Differentiable Loss Function for Time-Series[END_REF] is a smoothed formulation of DTW that computes the soft minimum of all alignment costs. Moreover, soft-DTW is a differentiable loss function, which makes it suitable for optimization tasks such as averaging time series and calculating centroids of time series clusters. The soft-DTW between two time series (a query sequence q and a reference sequence r), is defined in (3.15) considering the generalized min operator given in (3.16). Where ∆(q, r) is the cost matrix, A the alignment matrix, ⟨.⟩ the inner product and γ the smoothing parameter (γ ≥ 0). Note that the original DTW is obtained when setting γ = 0.

DT W γ (q, r) := γ min{⟨A, ∆(q, r)⟩} (3.15) γ min{a 1 , ..., a n } := min i≤n a i , γ = 0, -γlog n i=1 e -a i /γ , γ > 0.

(3.16)

Soft-DTW allows to estimate the barycenter or centroid rk of class C k , composed of a set of N k time series (r k,1 , r k,2 , ..., r k,N k ) of same length L, by minimizing the loss function defined in (3.17) through gradient descent.

min rk N k i=1 DT W γ (r k , r k,i ) (3.17)
Soft-DTW is not the only method in the literature that performs time series averaging. Another well-known method is DTW Barycenter Averaging (DBA), proposed by Petitjean et al. in [START_REF] Petitjean | A global averaging method for dynamic time warping, with applications to clustering[END_REF]. It is an iterative method that starts with an average sequence and improves it at each iteration following an expectationmaximization scheme. However, the best results in our case are obtained with soft-DTW, as shown in Fig. 3.34.

The figure presents the centroid estimation of signatures belonging to voltage sags caused by an upstream unbalanced fault (A2) using a set of reference signatures of the same class. We compare the centroid estimation results performed by (a) Euclidean distance, (b) DBA and (c) soft-DTW with γ = 0.5. Note that since the estimation is applied to the four dimensions simultaneously, the dimensions of the centroid signature are also time-correlated.

The results of the standard Euclidean averaging compared to the soft-DTW are significantly worse. The estimated centroid is obtained by averaging the curves without considering the temporal dimension. This is even more visible in the third component of the signature (voltage positive-sequence harmonic 2), where four peaks are observed instead of two peaks corresponding to the transients at the beginning and end of the sag. The soft-DTW also provides a smoother and more representative centroid than DBA: its characteristics are displayed in Table 3.5 for different values of gamma (γ). The best centroid estimation is obtained with soft-DTW with the lowest mean distance to the centroid for values of gamma γ ≤ 2. A correct setting of γ is required for an accurate centroid estimation. If γ is too low or too high, it can result in a distorted centroid, as observed in Fig. 3.35. An optimal value for γ should ensure a minimal distance between the reference signatures used for its calculation while maintaining a minimal distortion rate. For this, we calculate the distances D(r k,i , rk ) and the corresponding distortion rate Z (%) for each class in the synthetic data, as defined in equations (3.10) and (3.11) respectively. Fig. 3.36 illustrates the result, where the optimal value of γ opt = 0.1. The centroid estimation method can drastically reduce the total computational cost. However, it requires the correct setting of the parameter γ for an accurate centroid estimation. An incorrect value of this parameter can lead to non-representative and distorted class centroids, resulting in significant misclassification errors. Therefore, its applicability is conditioned to the choice of γ, for which an optimization is required.

Confidence score calculation

In addition to the class label of a new voltage sag, we propose to compute a confidence score associated with the prediction's result. Providing such a score increases the reliability of the classifier's output, as well as its interpretability. In the case of a decentralized approach to data processing, the availability of a confidence index could also reduce data traffic costs by limiting the amount of data uploaded to the server by selecting only those events with a low confidence index for further analysis. It can be noted that none of the methods in the literature described in Chapter 2 provide a confidence score associated with predicted label.

Confidence scores are commonly calculated for classifiers such as Support Vector Machines (SVM) and Artificial Neural Networks (ANN). These type of classifiers inherently provide a continuous output value due to their structure (softmax function for ANNs and distance to decision boundary for SVMs), which is used for calculating a posteriori membership probability. In addition, a confidence score calibration stage is usually recommended to provide probability estimates that are representative of the true likelihood, using post-processing techniques such as Platt's scaling [START_REF] Platt | Probabilistic Outputs for Support Vector Machines and Comparisons to Regularized Likelihood Methods[END_REF]. However, as discussed in subsection 3.6.2.1, k-NN and particularly 1-NN can not provide such scores. For this reason, we have extended the k-NN concept to a Nearest Neighborhood classifier presented in subsection 3.6.2.2 and we use the distance-to-class estimation for calculating a confidence score associated with the classifier's label prediction.

In this section we propose two distance-to-class-based confidence indexes: a probabilistic-based index (NB-KDE) and a relative distance-based index (RD). Both confidence scores are described hereafter, but their performance will be evaluated and compared in Chapter 4.

Probabilistic (NB-KDE) index

This first confidence index is obtained through a set of K binary Naive Bayes classifiers using a Kernel Density Estimator (NB-KDE). Each binary classifier is trained using a one-vs-rest approach. According to Bayes theorem in (3.18), the posterior probability P (y = C k |x) is proportional to the likelihood function estimated by the KDE P (x|y = C k ) and the prior P (C k ), with x = d k (q * ), and y being the label. Fig. 3.37 shows the likelihood function and the output score of a binary Naive Bayes classifier using KDE.

P (y = C k |x) ∝ P (x|y = C k )P (y = C k ) (3.18)
The probabilistic index NB-KDE associated with the prediction given in (3.14) is the output provided by the binary classifier k = y * , as defined in (3.19).

NB-KDE

= P (x|y = C y * )P (y = C y * ) (3.19)
The closer a new voltage sag signature q * is to the reference signatures r k of a given class C k in the database, the higher the NB-KDE confidence index will be. As expected, the database size influences the estimation of the likelihood function. An extensive and variate database will provide a reasonable estimate of the likelihood function by the KDE, thus, a more accurate NB-KDE confidence index. On the contrary, a reduced database will produce a poor likelihood function estimation, and the NB-KDE confidence index will not be considered reliable. The influence of the database size on the NB-KDE index will be presented in Chapter 4.

Relative distance-based (RD) index

The second confidence index is based on the calculated distances to the different classes. According to Ben-Israel et al. [START_REF] Ben-Israel | Probabilistic D-Clustering[END_REF], several relations can be assumed between the distance d k (q * ) and its membership probability p k (q * ), including the working principle defined in (3.20), where F (q * ) is a function depending only on q * . p k (q * )e d k (q * ) = F (q * ) (

This relation shows that the probabilities decrease exponentially as distances increase. Thus, the membership probability p k (q * ) can be defined as in equation (3.21), as proposed in [START_REF] Ben-Israel | Probabilistic D-Clustering[END_REF]. The confidence index RD is the membership probability p k (q * ) of the predicted class k = y * . A1). Although an extensive and diverse database would naturally improve the estimation of class distances d k (q * ), the relative distance-based index is in theory less sensitive to the database's size than the NB-KDE index. As for the NB-KDE index, the influence of the database's size on the RD index will be investigated in Chapter 4.

p k (q * ) = j̸ =k e d j (q * ) K i=1 j̸ =i e d j (q * ) , k = 1, 2, ..., K (3.21)

Conclusion

In this chapter, we presented a new methodology for classifying the causes of voltage sags based on multivariate time series analysis. It consists of four stages:

1. Data acquisition. Two data sources are available: numerical simulation (synthetic) data and real field measurements. For the generation of synthetic data, an industrial grid model is developed. Although the validation of the algorithm is to be performed on real field data, synthetic data is still necessary to cope with the limited size of the real datasets and to evaluate the algorithm's performance in different and controlled scenarios.

2. Pre-processing. This stage ensures the uniformity of the data in terms of length and sampling rate. It also allows for filtering incomplete voltage sags.

3. Feature extraction. The algorithm can transform voltage and current waveforms into electrically interpretable 4-dimension time series signatures through this process. This aspect is essential because it allows us to understand the overall decision-making process of the algorithm, which is a differentiating advantage compared to other methods in the literature. The description and the interpretation of the signatures belonging to the seven classes of voltage sag causes were also presented in detail.

4. Feature analysis. Using the DTW algorithm, we defined a distance measure for multivariate time series based on a space and time alignment process. Then, we proposed a distance-based algorithm for multivariate time series classification. It consists of an improved version of the well-known 1NN-DTW algorithm, whose major drawbacks are the sensitivity to outliers and high computational cost. Thus, we proposed a Nearest Neighborhood classifier, which instead of selecting the closest neighbor for the classification of a new signature, selects the closest class. The distance calculation of a new signature to a particular class can be obtained by estimating the mean distance value to all the class signatures or by estimating the distance to the classes' centroids. The latter significantly reduces the computational cost, but its applicability is conditioned to a correct parameter setting. Finally, we proposed to calculate two confidence indexes associated with the prediction.

In Chapter 4, we will compare the two methods for distance-to-class estimation and the two techniques for the confidence score calculation. We will also analyze the effectiveness of the feature extraction process by studying the separability between classes. The algorithm robustness will be evaluated in terms of noise levels and fundamental frequency variations, and the minimum database size will also be optimized. Finally, the generalization capabilities of the algorithm on synthetic and real field data will be presented.

Chapter 4

Classification of Voltage Sag Causes: Performance Analysis

Introduction

This chapter investigates the performance of the methodology presented in Chapter 3 in terms of class separability, sensitivity to noise, sensitivity to fundamental frequency variations, and computational cost. We will compare the two methods presented for the distance-to-class calculation: mean distance with bootstrapping and centroid estimation using soft-DTW. We will also compare the two methods for determining the confidence index associated with the prediction: NB-KDE index and RD index. Finally, the reference database size is optimized, and a cross-data source approach using synthetic and real field data is applied to evaluate the algorithm's global accuracy and generalization capabilities.

Class separability analysis

Evaluating the separability of classes is a way to verify the efficiency of the feature extraction process. We will study the separability of two classes by analyzing their intra-class and inter-class distances. The intra-class distance distribution H k is defined as the ensemble of distances D(x k,i , x k,j ) of all the pairs of elements x k,i and x k,j belonging to the same class C k . The inter-class distance distribution H k-k ′ corresponds to the ensemble of distances D(x k,i , x k ′ ,i ) between all the pairs of elements x k,i and x k ′ ,i from two different classes C k and C ′ k respectively. For this study, we use the Bhattacharyya coefficient (BC), as defined in equation (4.1), to evaluate the separability of the classes. This coefficient varies from 0 to 1, and can be interpreted as the overlap between two distributions H and H ′ .

BC(H, H

′ ) = x∈X H(x)H ′ (x) (4.1)
To study the class separability of a particular class α in a multi-class classification problem, we compute the Bhattacharyya coefficient BC(H α , H α-β ) between its intra-class distribution H α and the inter-class distributions H α-β , with β = {1, ...K}, α ̸ = β.

For example, in the case of two classes alpha and beta, they will be considered well separated if the overlap between the intra-class and inter-class distributions is close to zero. This would correspond to BC(H α , H α-β ) and BC(H β , H β-α ) close to zero. For all the considered seven classes, the BC coefficients are computed and displayed in Table 4.1.

From these results, we can consider that the classes are sufficiently well separated: all the BC values are close to zero, except on the diagonal, where they are equal to one because they correspond to the calculation between two identical distributions. The largest overlap (worst case) corresponds to BC = 0.07, for the intra-class distribution H C2 and the inter-class distribution H C2-C1 . Indeed, the first three ISCs of their signatures are identical. Only the current component (fourth ISC) discriminates the upstream event from the downstream one. The overlap between these two distributions is due to the higher dispersion of

Table 4.1: Bhattacharyya coefficient (BC) between the intra-class and inter-class distance distributions

Inter-class distribution H α-β H α distribution Intraclass H α-A1 H α-A2 H α-B1 H α-B2 H α-C1 H α-C2 H α-D H A1 1
H C2 compared to H C2-C1 .

Evaluation of the performance using synthetic data

This section compares the performance of the two distance-to-class estimation methods, and the two confidence indexes presented in Chapter 3. We analyze the sensitivity of the different methods to increasing noise levels and variations of the fundamental frequency. We will study how to determine the optimal size of the database and compare the computational costs. For this analysis, we use the synthetic dataset presented in the previous chapter to evaluate the algorithm's limits and validate the methodology's relevance.

Sensitivity analysis

Metrics and criteria

We can evaluate the sensitivity to noise and frequency variations by analyzing their impact on the classification's accuracy, and the confidence in the prediction. The standard metrics for classification algorithms commonly used to evaluate the performance of label prediction are precision, recall, and F1-score2 . The F1-score, described in equation (4.2), is a good metric for summarizing the first two, but it is mainly applied to balanced datasets. A perfect classification is obtained with an F1-score equal to one.

F 1 = 2 * precision * recall precision + recall (4.2)
However, the metrics mentioned above are not suitable for evaluating the confidence score associated with the prediction provided by the algorithm. A more appropriate metric for this task is Log-Loss. This metric penalizes predictions with a low confidence score. The definition of Log-Loss is given in equation (4.3), where y = {0, 1} and p is the associated probability estimate with p = P (y = 1).

L log = -(ylog(p) + (1 -y)log(1 -p)) (4.3)
Although raw Log-Loss values can be hard to interpret, lower values mean predictions with higher confidence. For instance, a perfect classifier would have a Log-Loss equal to 0. A random guess Log-Loss baseline score can be useful for interpreting this metric. The Log-Loss value for p = 0.5 is L log = 0.693. Any value higher than this baseline (represented as a red dashed line on figures 4.3 to 4.7) can be interpreted as worse than random guessing.

Sensitivity to noise

To study the classifier's robustness to noise, we constitute two different reference databases used for "training"3 , and the algorithm is tested on 16 test sets with increasing noise levels.

• Reference database A contains the original 700 synthetic voltage sags from the synthetic dataset (with no added noise)

• Reference database B comprises the original 700 synthetic voltage sags plus 700 sags with additional white Gaussian noise at SNR=25 dB (1400 sags in total) Reference database B is used to evaluate if enriching the "training" set with noisy data can improve the classification results. It should be noted that reference database B is enriched with noisy signatures at SNR = 25 dB, which is a noise level not included in any test set. The closest are test set no. 8 and no. 9 with SNR = 26 dB and SNR = 24 dB, respectively.

The results of the experiments with reference databases A and B, for the distance-to-class calculation method using bootstrapping, are illustrated in Fig. 4.3, and for the centroid estimation method in Fig. 4.4.

The evolution of the F1-score (Fig. 4.3a and Fig. 4.4a) is almost identical for both methods, independently of the used reference database. F1-score reaches its maximum value close to 100% from 40 dB to 20 dB, and slowly decreases in the range 20 to 10 dB to reach 90%. The main differences are found in the estimated confidence indexes. Independently of the method or the reference database, the mean values of the NB-KDE ). However, for SNR levels between 30 to 10 dB, we note a rapid degradation of the NB-KDE index, which exceeds the random-guess baseline around 20 dB for the bootstrapping method in Fig. 4.3b, and around 25 dB for the centroid estimation method in Fig. 4.4b. The use of an enriched database (reference database B) helps to improve the results of the NB-KDE index. This is an expected behavior since a more diverse database leads to a larger likelihood distribution. However, the performance of this index within this range of noise levels is still very poor.

In contrast, the RD index seems much more stable, with an average value above 75% for the bootstrapping method in Fig. 4.3a, and up to 80% for the centroid estimation method in Fig. 4.4a, even for noise levels corresponding to SNR of 10 dB. We also note that the addition of noisy data to the reference database (reference database B) does not significantly influence the RD index. In summary, we can conclude that the classification's performance represented by the F1-score for both distance-to-class methods is robust enough in the standard range of noise levels corresponding to SNR between 40 to 15 dB (above 95%). This is because the feature extraction step only retains the low-frequency components, implicitly filtering out the high-frequency noise when applying STFT. Fig. 4.5 illustrates the second ISC (voltage negative-sequence harmonic 1) of a sag caused by an upstream balanced fault with different noise levels. It can be observed that the noise in the voltage waveforms have a low impact on the signatures. Finally, we also note that the RD-index is significantly more robust than the NB-KDE index for noise levels with an SNR between 25 and 10 dB.

Sensitivity to fundamental frequency variations

The nominal supply frequency in France is 50 Hz. However, this frequency can slightly fluctuate around this value due to changes in the supply and demand for electricity in the grid. To evaluate the sensitivity to frequency variations around the fundamental frequency (50 Hz), we proceed as before with experiments based on two reference databases and 6 test sets with increasing levels of frequency variation.

• Reference database A contains the original 700 synthetic voltage sags of the synthetic dataset at the rated fundamental frequency of 50.0 Hz

• Reference database C contains the 700 original synthetic voltage sags plus 105 sags at 49.75 Hz and 105 sags at 50.25 Hz (910 sags in total)

• Test sets, 6 in total. Each test set contains 210 synthetic voltage sags with a fundamental frequency F = 50Hz ± ϵ, where ϵ = {0.0, 0.

The range of frequency variation corresponds to the maximum frequency fluctuations (50Hz ±1%) allowed by the French and European regulatory standards regarding the power supply at the distribution level for synchronous connection to an interconnected system [START_REF]for Electrotechnical Standarization[END_REF].

The results for the bootstrapping distance estimation method are plotted in Fig. 4.6, and for the centroid estimation method in Fig. 4.7.

For reference databases A and C, the values of the F1-score are above 95% for the bootstrapping method in Fig. 4.6a and the centroid estimation method in Fig. 4.7a.

The NB-KDE index is more sensitive than the RD index to increasing frequency variation levels for both methods. However, its mean and Log-Loss values are improved using the enriched reference database C. This improvement is more significant with the centroid estimation method illustrated in Fig. 4.7a, where the mean value at ϵ = 0.2 changes from 50% with reference database A, to almost 90% with reference database C (+40% improvement). An improvement in the bootstrapping method (+20% for ϵ between 0.2 to 0.5) is also visible in Fig. 4.6a. The Log-Loss values for this index for both methods evolve similarly (Fig. 4.6b and Fig. 4.7b). The use of reference database C helps to noticeably decrease the Log-Loss values for the bootstrapping and centroid estimation methods, with values close to the random guess baseline for ϵ = 0.3. However, the overall results obtained with the NB-KDE index are not sufficient. On the other hand, we note that the RD index is significantly more stable, remaining above 90% even for fluctuations up to 0.5 Hz, for both distance-to-class methods and both reference databases (Fig. 4.6a and Fig. 4.7a). Regarding the Log-Loss values in Fig. 4.6b and Fig. 4.7b, the RD index remains low and stable despite frequency variations, with values significantly lower than the random guess baseline. We note that the enriched reference database C has a more noticeable impact on this index in the case of the centroid estimation method (Fig. 4.7). Reference database C slightly improves the RD index's results in terms of mean and Log-Loss values.

The impact of fundamental frequency variations on the signatures is illustrated in Fig. 4.8, with the second ISC (voltage negative-sequence harmonic 1) of a sag caused by a transformer energizing. We observe that although the frequency variation causes visible oscillations in the signature, the global shape that allows the classification is still preserved. Finally, we can conclude that the classification performance measured by the F1-score is very high (above 97%) even for fundamental frequency variations of ϵ = 0.5Hz. We also note that the RD index proves to be more robust once again, achieving a mean value higher than 90% for both distance-to-class methods and both reference databases.

Minimum database size setting

Classification algorithms in the literature generally require large amounts of data to perform, from a few hundred to thousands of training samples. Unfortunately, in industrial applications, the availability of disturbed electrical data in various operating conditions is hard to obtain. Therefore, in this area of research, it is essential to perform the classification with a small database size. This section determines how to adjust the optimal size of the database to achieve the best possible compromise in terms of performance.

The results are obtained using the original synthetic dataset of 700 sags (no added noise or fundamental frequency variations). The training database contains a determined and increasing percentage of the synthetic dataset, and the test set contains the remaining data. We also use a random permutation strategy of five balanced splits for more representative results.

Results for the mean distance estimation with bootstrapping are illustrated in Fig. 4.9, and for the centroid estimation method in Fig. 4.10.

A good trade-off between minimum database size and high confidence for the NB-KDE and RD indexes is obtained for a database size containing no more than 20% of the original synthetic dataset. The F1-score and the RD index reach high scores for even less data, but a minimum of 20 signatures per class are required in the database to obtain a reliable score for the NB-KDE index. Above this ratio, there are no significant improvements in performance. We present hereafter the detailed classification results for a database with 20 signatures per class (140 signatures for the 7 defined classes). The reference database contains 20% of the original synthetic dataset and the remaining 80% is used for testing. The results for the mean distance with the bootstrapping method are presented in Table 4.2 and for the centroid estimation method in Table 4.3. They are obtained through a cross-validation strategy using five balanced splits, so the whole synthetic dataset is evaluated. The results for the label prediction, reflected by the F1-score, are optimal for both distance-to-class estimation methods with a 100% success rate for each class. In terms of confidence in the prediction, the best results are achieved with the bootstrapping method combined with the NB-KDE index (99.35% of mean confidence across all classes). However, the sensitivity analysis in Section 4.3.1 showed that the NB-KDE index is highly sensitive to increasing noise levels or frequency variations. On the other hand, the RD index is slightly lower but still with very satisfactory results for both distance estimation methods: 96.27% of mean confidence for the bootstrapping method and 96.46% for the centroid estimation method. Therefore, we select the RD index as confidence score in the following sections.

Computational cost evaluation

Using the synthetic dataset, we compare the computational cost of the distanceto-class calculation methods: mean distance estimation with bootstrapping and centroid estimation using soft-DTW. The execution time for classifying a single signature (test runtime) with different database sizes is presented in Table 4.4. In addition, the computation time used by the centroid estimation method to estimate the barycenters is also given as a reference (training runtime). The bootstrapping method does not require any prior calculations before implementing the algorithm, therefore its train runtime is not provided. Note that only the label prediction time is measured to compare the speed performance of these two methods. The time to compute the confidence index is not considered, as it is negligible compared to the time to compute the distance-to-class. The experiments are performed with an Intel Core i7-8750H processor in a Windows 10, 64 bits operating system. Ten successive runs are performed and the average value is provided. The centroid estimation method is faster than the mean distance estimation with bootstrapping approach during the test runtime. The test time ratio of this method is proportional to the number of events per class N k . By reducing the reference signature database to a single centroid per class, we can efficiently reduce the overall classification computation time by a factor equal to N k . The centroid estimation testing time is stable whatever the number of signatures per class N k . The centroid estimation training runtime increases linearly with the size of the database. However, the centroid estimation is only performed once, prior to the implementation of the algorithm.

In the previous sections, we studied the performance of two methods for distanceto-class estimation and two methods for calculating the confidence score associated with the prediction. The sensitivity analysis to noise and frequency variations highlighted the robustness of the RD index compared to the NB-KDE index. For a fixed-size database, the mean distance estimation with the bootstrapping method provided similar results to the centroid estimation with soft-DTW. However, the latter outperforms the first in terms of computational speed by a factor equal to N k (number of signatures per class). The centroid estimation with soft-DTW and RD index as confidence score are the most suitable methods in terms of speed and robustness. Therefore, they are retained for validation using real field data in the following sections.

Validation using real field data

The performance of our proposal with field data is evaluated in this section. We also evaluate its capability of generalization with different data sources (synthetic and real field data). The analysis is performed under the scope of the confusion matrix 4 and F1-score metrics. Misclassification errors are also assessed, and the results highlight the usefulness of the proposed confidence index in identifying and excluding possible false results.

Cross-data source evaluation

The overall classification accuracy of some methods in the literature is very high, (over 99%) [START_REF] Li | Adaptive S transform for feature extraction in voltage sags[END_REF][START_REF] Yikun | Classification Method of Voltage Sag Sources Based on Sequential Trajectory Feature Learning Algorithm[END_REF][START_REF] Zheng | Recognition method of voltage sag causes based on two-dimensional transform and deep learning hybrid model[END_REF][START_REF] Zheng | Recognition method of voltage sag causes based on Bi-LSTM[END_REF][START_REF] Wang | Method of Voltage Sag Causes Based on Bidirectional LSTM and Attention Mechanism[END_REF]. Nevertheless, there is an important limitation regarding the required data for the development and deployment of these algorithms: these classifiers are systematically trained and tested using data from the same data pool (either simulation or real data).

In the case of deep learning classifiers, a mix of synthetic and real data is used for training to deploy the algorithm on real data [START_REF] Zheng | Recognition method of voltage sag causes based on two-dimensional transform and deep learning hybrid model[END_REF][START_REF] Wang | Method of Voltage Sag Causes Based on Bidirectional LSTM and Attention Mechanism[END_REF][START_REF] Zheng | Recognition method of voltage sag causes based on Bi-LSTM[END_REF]. This is partly because these algorithms require important amounts of data to converge, and this strategy is used to accelerate the learning process and improve overall results. However, the authors have not reported the results of these classifiers when exclusively trained with simulation data and then tested on real data. We know from [START_REF] Bollen | Classification of Underlying Causes of Power Quality Disturbances: Deterministic versus Statistical Methods[END_REF] and [START_REF] Axelberg | Support Vector Machine for Classification of Voltage Disturbances[END_REF] that this is not a trivial task. Both authors report good results when training an SVM classifier with field data from one network and testing it on a different network. However, the classifier's performance significantly deteriorates when trained with simulation data and tested with field data. The authors also pointed out that the usefulness of these methods for commercial purposes depends on the ability of the classifier to be factory pre-trained and perform accurately when deployed in different networks, as it is not realistic for customers to train the algorithms themselves. Furthermore, from an implementation perspective, it is a strong demand for classifiers to be based solely on synthetic data, as it is more accessible than field data [START_REF] Axelberg | Support Vector Machine for Classification of Voltage Disturbances[END_REF].

We evaluate the accuracy and generalization capabilities of the classifier with the centroid estimation method. For this, we conducted four experiments described in Table 4.5, using the synthetic and real datasets described in Chapter 2. The synthetic and real datasets are split into two groups: 20 signatures per class for the reference database (training), and the rest for testing.

Note that for the real dataset from three industrial sites, classes B1 and C2 are not present since no recordings were available during the monitoring period of time (see Table 3.2). Also note that class C1 has very few signatures: the real field reference database comprises 5 events instead of 20 for the other classes. For experiment 2, the results in Table 4.7 show that the F1-score is equal to 90.19% for 5 classes out of 7, and the average confidence score is 92.64%. The errors are due to the misclassification of 4 voltage sags belonging to class A2 (upstream unbalanced fault) predicted as events of class C1 (upstream transformer energizing).

The classification errors in experiment 2 are due to a distorted estimated centroid for class C1, which caused a forced time alignment between four events of class A2 with the centroid of class C1. The distorted centroid for class C1 is mainly due to one out of five signatures used for the estimation that presents some oscillations and irregularities. Since the number of signatures used for the centroid estimation is low (only five), a single outlier can affect the barycenter estimation process. This problem could be prevented by: discarding outliers from the reference signatures used for centroid estimation, by increasing the number of reference signatures, or by increasing the value of γ to obtain a smoother centroid.

The most significant results are those from experiments 3 and 4, since the reference training database and the test set belong to different data sources.

The results of experiment 3 are displayed in Table 4.8. They show that the algorithm can use a synthetic database for training and obtain satisfactory results when tested on real field data with a 99.32% average F1-score and an average confidence score of 87.02% for 5 out of 7 classes. The results from this experiment are the most relevant since, as previously mentioned, an algorithm with commercial purposes should be ideally trained only with synthetic data before deploying it in actual grids to evaluate real field data. In experiment 3, only two sags of class A2 (upstream unbalanced line fault) were misclassified as belonging to class A1 (upstream balanced line fault). This classification error is due to an overlap of the two peaks in the third ISC of the signature (voltage positive-sequence harmonic 2) because of the short duration of both sags (less than 30 ms). This overlap is responsible for a poor time alignment with both reference signatures A1 and A2. Fig. 4.11 and Fig. 4.12 show a comparison between the new sag signature (30 ms) and a reference signature of class A1 (50 ms), and a reference signature of class A2 (50 ms), respectively.

We can observe that an incorrect time misalignment led to the distortion of the signature in both cases, resulting in a distance-to-class A1 of D = 5.52e -02 (false class) and D = 5.96e -02 to class A2 (true class). However, both distances are significantly higher compared to the average inter-class distances for class A1 and A2 (lower than D = 3.50e -02). This is a limitation of the feature extraction process for sags with a duration inferior to 30 ms. However, classification errors due to poor time alignment can be detected and may raise an alert if the confidence score (RD index) is low and the intervention of a human expert is required. In this case, a human eye would have rapidly observed that the signature of the new event is closer to the signature of class A2 than to class A1, due to the voltage negative-seq H1 (second ISC), despite the peak overlap in the third ISC. Error analysis using the confidence index will be further studied in the next section. Finally, although experiment 4 is not to be replicated in practical applications, the results presented in Table 4.9 are still interesting from a validation perspective for the proposal's capability of generalization. For this experiment, the algorithm achieves a F1-score of 100% and 93.23% for the RD index.

These last two experiments demonstrate the algorithm's capacity to properly generalize across data from different sources (either synthetic or real). It should be noted that such tests and have not been performed by other methodologies in the literature for the classification of voltage sag causes.

Accuracy and error analysis

This section analyzes the classification errors for experiment 3 (synthetic database and real test set) presented in 4.8. Two voltage sags of the real dataset belonging to class A2 were misclassified as belonging to class A1. If we analyze the relative distance-based confidence index of all the voltage sags, we observe that only 5 events out of 385 were classified with a confidence index lower than 60%, as illustrated in Fig. 4.13. Among these events, two correspond to misclassification errors.

The classification error appears between two classes that only differ on the balanced/unbalanced nature of the event. Classes A1 and A2 correspond to voltage sags due to upstream line faults. The consequences of this error could even be considered minor since the location (upstream) and the source event (line fault) of the sag are the most helpful information for the industrial client. A prediction with a confidence index below a certain threshold would trigger an alert for further analysis by an expert. The end-user of the system can set this threshold according to his needs. A threshold set at 60% for the RD index seems to be a good compromise between the number of events triggering an alert and the risk of overlooking a classification error. The possibility of setting such a threshold reduces the time spent by experts in analyzing voltage dips. For the evaluated data, only 5 events out of 385 (i.e. less than 2%), would have required a detailed study.

The two events classified in class A1 have an RD index of 41.01% and 42.34%, respectively. These values are significantly lower than the median RD index for this class, close to 91%. In addition to this, we note that the next highest membership probability for both events are 29.66% and 38.93%, respectively, which correspond to class A2, the true class. We can also observe that some classes with low membership probability can be quickly excluded. In the case of a more detailed analysis required by an alert, this information can help reduce the number of possible class labels to be analyzed from seven to only two or three.

Conclusion

This chapter analyzed the algorithm's performance in terms of class separability, robustness to noise, and fundamental frequency variations for classifying the causes of voltage sags. Two methods for distance-to-class estimation were compared, and the centroid estimation with soft-DTW method was found to be as efficient but significantly faster than the bootstrapped mean distance estimation method. Similarly, two methods for calculating the prediction confidence index were investigated, and the RD index was found to be much less sensitivity to noise and frequency variations than the NB-KDE index. Table 4.10 briefly summarizes the criteria for the selection. The results proved that the proposal is resilient for noise levels up to SNR = 15 dB and fundamental frequency variations up to ϵ = ±0.5 Hz. The algorithm was evaluated using a cross-data source approach using synthetic and real field data through four experiments, showing very good generalization capabilities. The algorithm reached a F1-score of 100% across seven classes with a reduced synthetic database of only 140 events when tested on synthetic data. However, the most relevant results were obtained with a synthetic signature database for training and testing on real field data collected from three different industrial sites. The algorithm reached a F1-score of 99.32% for five out of the seven defined classes.

In summary, the main advantages of the proposed algorithm are:

1. The reduced amount of data necessary to build the signature database 2. The signature reference database can be entirely composed of synthetic data 3. The good generalization capabilities when implemented on real field data, even for different industrial sites 4. The electrical intepretability of the signatures and the decision-making process

The computation of a confidence index associated with the prediction

The first three characteristics make the algorithm easy to implement in real industrial applications with no previously recorded data. The last two characteristics make the troubleshooting process easier and increase the general interpretability of the decision-making process, which is a demand from industrial customers for reliability issues.

Chapter 5

Impact of Voltage Sags in Industrial Grids

Introduction

Once a voltage sag has been detected by the monitoring device, and its root cause has been identified, the next step is to analyze the impact on industrial equipment.

Voltage sags are indeed one of the main causes of load self-disconnection [START_REF] Yamashita | Identification of Individual Load Self-disconnection Following a Voltage Sag[END_REF]. The impact in terms of load-shedding caused by a voltage sag depends on the characteristics of the sag and the resilience of the industrial grid and its processes. This study aims to propose a methodology capable of estimating the self-disconnected load composition in an industrial facility after a voltage sag. For this, we consider that we only have access to the electrical measurements provided by a unique monitoring device at the LV side of the industrial grid's main MV/LV transformer.

In this chapter, we present the problem of load self-disconnection after a voltage sag and introduce some of the most common methods in the literature for load estimation, particularly for load composition estimation. Then, we propose a first approach for the problem of self-disconnected load composition estimation. We present the obtained results for a simple case study and analyze the method's limitations. Finally, we discuss the challenges and perspectives.

Self-disconnected loads due to voltage sags

Voltage sags may or may not cause an impact on the industrial facility. A voltage sag with impact refers to a sag causing load self-disconnection of sensitive equipment. The duration of the impact is variable and depends on the time needed to restore the normal operation of the affected processes. Indeed, the restoration can be automatic with certain variable speed drives restarting the motors once the voltage conditions are met, or may require human intervention for more complex processes. The total duration before the restoration of the stopped processes has a direct impact on the financial losses, as presented in Chapter 1. .1a presents a voltage sag with a magnitude of 87%, without impact. We can observe that the current and the active power return to their pre-sag values after the voltage is restored. In the contrary, Fig. 5.1b illustrates a voltage sag with a magnitude of 43%, causing the disconnection of multiple loads. The impact is particularly visible when comparing the pre-sag and post-sag values in current (∆I) and active power (∆P ). This voltage drop causes a loss of 64% of the total active power consumed by the site.

This study aims to develop an algorithm that provides as much information as possible on the disconnected loads following a voltage sag, with limited access to electrical measurements. However, it is not feasible to determine the actual characteristics of the devices connected downstream to the single monitoring point. Therefore, we group the loads based on their dynamic response in five load categories including the most common equipment in industry [START_REF] Milanovic | Methodology for Estimation of Dynamic Response of Demand Using Limited Data[END_REF]:

1. Direct connected induction motors 2. Adjustable speed drivers (ASD) 3. Resistive loads 4. Switched-mode power supply (SMPS) feeding electronic devices

Energy efficient lighting

The objective is to estimate within these five load categories the composition of the disconnected loads due to a voltage sag.

Load composition estimation is actually an area of application of load estimation techniques. So, in the following, we first introduce a brief literature review on load estimation techniques. Then, we propose a first approach to estimate the self-disconnected load composition after a voltage sag.

Load estimation methods: Literature review

Accurate load estimation is essential for power analysis, planning and control [START_REF] Arif | Load Modeling-A Review[END_REF]. Load estimation can be divided in two stages: load modeling and parameter estimation.

Load models

Load models can be designed with analytical relations (polynomial, exponential, non-linear, etc.) that describe the load's response to variations of the input voltage or/and frequency. Models can also be built from equivalent electrical circuits, whose parameters are related to the characteristics of the load.

Machine learning techniques like Artificial Neural Networks (ANN) can also be used to model loads. They are helpful for modeling loads whose physical structure is unknown or whose electrical behavior is too complex to be formulated through mathematical equations. However, as discussed in Chapter 2, these methods require huge amounts of data for training. Their accuracy and particularly their generalization capabilities will depend largely on the diversity and representativeness of the data used for training.

The domain of validity should be considered when selecting a load model [START_REF] Renmu | Composite Load Modeling via Measurement Approach[END_REF]. Sometimes the model is an approximation around an operating point, for instance up to 5% voltage sag [START_REF] Milanovic | Methodology for Estimation of Dynamic Response of Demand Using Limited Data[END_REF]. However, in case of a deeper voltage sag, several electrical components may exhibit a very different behavior than the one observed in nominal conditions, and in consequence, the model becomes invalid.

The final use of the model is also a criterion that should be considered [START_REF] Renmu | Composite Load Modeling via Measurement Approach[END_REF]. For example, stability analysis, electrical consumption forecasting, load parameter estimation, load composition, etc. For example, a steady-state model of aggregated loads is sufficient for consumption forecasting studies. However, dynamic and physics-based models are required to estimate the electrical and mechanical parameters of electrical motors.

The selection of a load model depends on the characteristics of the parameters to estimate, which can be:

• The coefficients of analytical equations, with no direct physical interpretation • The electrical or mechanical parameters of a specific device • The load composition1 of a residential, industrial or commercial grid In summary, load models can be grouped into three categories: static models, dynamic models and composite models.

Static models

Static models represent loads that reach a new equilibrium point immediately after a voltage disturbance without significant transient. They are suitable for analyzing steady-state characteristics. The most common models in this category are [START_REF] Arif | Load Modeling-A Review[END_REF]:

1. ZIP model: it consists of two polynomial equations describing the active and reactive power as defined in (5.1). Each equation is represented by three components, namely, a constant impedance part (Z), a constant current part (I), and a constant power (P) part. P 0 and Q 0 are the active and reactive 2. Complex Load Model (CLOD): is a model that includes: large and small motors, a discharge lighting block, a transformer saturation block, a constant power load (MVA), shunt capacitors, and a series impedance. It has been included in the Siemens PTI PSS/E stability program and successfully implemented for modeling the dynamic response of a pulp and paper mill facility to a voltage sag in [START_REF] Li | Dynamic load modeling for industrial facilities using template and PSS/E composite load model structure CLOD[END_REF], after conducting an in-depth load survey for parameter setting. [START_REF] Arif | Load Modeling-A Review[END_REF]. In [START_REF] Liang | A New Composite Load Model Structure for Industrial Facilities[END_REF], the author proposed a generic composite load model structure for industrial facilities using a similar structure as the WECC composite model. The method falls into the component-based approach category, which means that the parameters of the model are defined based on survey results.

In summary, static models are relevant for steady-state studies such as the load consumption estimation of aggregated loads in residential or commercial areas. However, they are not adapted for stability analysis as they ignore the dynamic behavior of the loads, whereas dynamic models are preferred. Loads can also be modeled using analytical equations or physics-based models. Analytical equations are useful if we are only interested in replicating the load's behavior. However, these models are not adapted if the goal is to model a group of loads as these equations do not allow to differentiate them.

Estimation methods

There are two main approaches for estimating load model parameters: componentbased and measurement-based methods [START_REF] Arif | Load Modeling-A Review[END_REF].

Component-based methods

Component-based methods rely on general knowledge of physical behavior and mathematical relationships describing loads. This approach is mainly used for power flow, voltage profile, and load demand studies [START_REF] Collin | Component-based aggregate load models for combined power flow and harmonic analysis[END_REF]. These types of studies do not require highly precise models as long as the global response of the model is representative of the real installed load. The load model parameters are set according to tables with standard values [START_REF] Liang | A New Composite Load Model Structure for Industrial Facilities[END_REF]. Nonetheless, such information is not always available or precise enough.

Measurement-based methods

Measurement-based methods are used to develop accurate load models. They use identification and estimation techniques to fit the structure and parameters of the model with measurement data. Some of the most usual algorithms for parameter identification are: least-squares [START_REF] Liu | The load modeling and parameters identification for voltage stability analysis[END_REF], genetic algorithms [START_REF] Bai | A Novel Parameter Identification Approach via Hybrid Learning for Aggregate Load Modeling[END_REF], fuzzy regression [START_REF]Advanced Load Modeling[END_REF][START_REF] Sagi | Evaluation of a load composition estimation method using synthetic data[END_REF] and particle swarm optimization (PSO) [START_REF] Regulski | Estimation of Composite Load Model Parameters Using an Improved Particle Swarm Optimization Method[END_REF].

In our application, field measurements provided by a monitoring device are available at the LV industrial site. Therefore, the use of measurement-based methods are relevant to estimate the load composition. However, in the load estimation domain, very few methods are focused on load composition estimation. Most of the measurement-based approaches are used for electrical/mechanical parameter estimation of specific devices and equivalent load response modeling for stability analysis [START_REF] Renmu | Composite Load Modeling via Measurement Approach[END_REF][START_REF] Regulski | Estimation of Composite Load Model Parameters Using an Improved Particle Swarm Optimization Method[END_REF].

The closest study to our problem is reported in a technical report published by the Electric Power Research Institute (EPRI) on a methodology for load composition estimation [START_REF]Advanced Load Modeling[END_REF]. The study is based on a composite model composed of: three ASM with different power ratings, an incandescent and discharge lighting, a thermostatic load and an electronic load. The active and reactive powers during a voltage sag are used as input of a fuzzy regression algorithm. The authors achieved good approximation with simulation data (a mean estimation error less than 2%). However, the models of individual loads are very simplified analytical equations that do not accurately reproduce their behavior in real conditions. For validation, only qualitative results were provided with field data. Sagi et al. [START_REF] Sagi | Evaluation of a load composition estimation method using synthetic data[END_REF] proposed a method for load composition estimation using fuzzy regression and evaluated the proposal with synthetic data. The authors used a component and physics-based model that included: incandescent and discharge lighting, electronic devices, ASD, thermostatic, small, medium and large motors. The authors also used generic load models assuming standard parameter values. The results showed that load composition could be estimated using the active and reactive power information with an acceptable accuracy despite a parameter mismatch (particularly power factor mismatch). The estimates for individual motor types (small, medium and large) were mismatched by a maximum of 8%, but the estimate for motor loads as a group was more accurate (less than 2% estimation error). Soon et al. [START_REF] Lee | A Reduced Multivariate-Polynomial Model for Estimation of Electric Load Composition[END_REF] proposed a methodology based on the harmonic content of individual loads and a reduced multivariate-polynomial model for the estimation of electric load composition. However, the instantaneous harmonic current signatures i(t) of the four considered loads (incandescent lighting, fluorescent lighting, computer, and motor drive) are modeled using static analytical equations independent of voltage variations.

Finally, Duan et al. [START_REF] Duan | Neural network approach for estimation of load composition[END_REF] proposed a ANN to estimate the load composition of commercial and residential grids. The algorithm exploits the current harmonic content. It is based on the assumption that the total measured current waveform can be approximated as the weighted sum of a set of current waveforms flowing in the connected loads. However, the harmonic content of each load is considered as constant and is not obtained from a load model. The authors use as input the harmonic spectrum at the service point and a table of the harmonic content of typical residential electrical loads (9 in total). The output are the weights of the individual loads (load composition). Only two experimental results were provided, with a percentage relative error (PRE) of 9% in both cases.

The load composition methods found in the literature either use the active and reactive powers, or the current and its harmonic content. The methods using active and reactive powers use load model approximations that are only valid for shallow voltage variations (around ± 5% of nominal value). These models are not accurate enough to replicate the response of the electrical devices to voltage sags (between 90% and 10% of the remaining voltage). In the other hand, methods using the current and its harmonic component provide new and useful information for the load composition estimation, particularly for differentiating devices with power electronics. However, the harmonic current models are based on static analytical equations that do not take into account voltage variations.

Self-disconnected load composition estimation:

A first approach

Inspired by the load estimation methods in the literature, we propose an approach based on two steps to estimate the composition of self-disconnected loads due to voltage sags:

Load composition estimation after a voltage sag without impact

We first consider voltage sags without impact (no self-disconnected loads).

• We have access to a detailed list of the electrical equipment in the industrial facilities and their characteristics, and therefore it is possible to model each one of the existing devices. However, their operating cycles are unknown.

• The voltage tolerance curves (VTC) are known and are similar for all the devices in the same category. In addition, the VTCs have a rectangular shape with two parameters: V max and t max .

• The voltage sags are caused by upstream line faults (balanced or unbalanced) in a parallel line feeder.

Load modeling

To estimate the load composition of an industrial grid, we should have for each load category an accurate model that is able to reproduce the behavior of the group, before, during and after the voltage sag. Thus, we must use dynamic models.

The use of models based on analytical equations requires setting the coefficients beforehand. However, our goal is to estimate the weights corresponding to the load composition. Therefore, dynamic analytical models, which can be very complex, are not the best choice for our problem. Consequently, we select the simulation software EMTP-RV to have access to more accurate physics-based models. The voltage and current measurements during the sag are recorded by the monitoring device in the actual industrial grid. Then, the three-phase voltages are injected into the load models of EMTP-RV. Finally, the currents absorbed by each load during voltage sags are collected and used by the load composition estimation algorithm.

Load composition estimation

The proposed load composition estimation method is a measurement-based approach. The constrained least-squares algorithm is used to estimate, at the voltage sag occurrence, the weighting coefficients that represent the load composition, at first without load tripping.

The dynamic response of each load category to a voltage variation, such as a voltage sag, is different. The first two equations of the load models should express the relation between the total active (P) or reactive (Q) power variations during a sag as the sum of the powers drawn by each load.

The harmonic content of power electronic devices such as ASD, three-phase rectifiers, and energy-efficient lighting is an important characteristic. One may intuitively think of the distortion power (D) as the third equation for load estimation. However, it is important to note that harmonics vary both in amplitude and phase. Therefore, the total distortion power (D) at the main feeder is not equal to the algebraic sum of the distortion powers drawn by each load. Therefore, harmonics can only be added in their complex form. We retain the fifth harmonic current as a reference for the harmonic content. This harmonic is usually the second component with the highest amplitude after the fundamental due to the presence of six-pulse rectifiers in the loads with power electronics.

The load composition is obtained by the estimation of the weighting coefficients that represent the different load categories. The equations are described in (5.5) and (5.6), where the error estimate e = {e P , e Q , e I } is minimized. P X (t), Q X (t) and the real component of Ih5 X (t) are the measurements per phase of length L, with X = {A, B, C} obtained from by the monitoring device in the actual industrial grid. PX (t), QX (t) and Î h5 X (t) are the load model's estimates. They are calculated by adding the individual response of each load defined as Pi,X (t), Qi,X (t) and Î h5 i,X (t), multiplied by a weight coefficient α i,X . To ensure the proportionality between the active power, reactive power and fifth harmonic current, the equations are normalized with the nominal active power prior to the sag P X (0), PX (0) and Pi,X (0).

                       e P = L t=0 [P X (t) -PX (t)] 2 e Q = L t=0 [Q X (t) -QX (t)] 2 e I = L t=0 [Ih5 X (t) -Î h5 X (t)] 2 (5.5)                        PX (t) PX (0) = 5 i=0 α i,X Pi,X (t) Pi,X (0) QX (t) PX (0) = 5 i=0 α i,X Qi,X (t) Pi,X (0) Î h5 X (t) PX (0) = 5 i=0 α i,X Î h5 i,X (t) Pi,X (0)
(5.6) Fig. 5.4 illustrates the measurement-based approach, where the objective is to determine the weight coefficients that minimize the error e between the measurements of the actual industrial grid y and the model estimates ŷ.

Figure 5.4: Measurement-based approach for load composition estimation

There are fifteen weighting coefficients to estimate for the five load categories in three phases. However, except for lighting, all of the modeled loads are threephase loads. Thus, for these loads we can consider that: α i,A = α i,B = α i,C = α i , which reduces the number of variables to be estimated to seven, as described in (5.7).

α =       α 1,A α 1,B α 1,C α 2,A α 2,B α 2,C α 3,A α 3,B α 3,C α 4,A α 4,B α 4,C α 5,A α 5,B α 5,C       =       α 1 α 1 α 1 α 2 α 2 α 2 α 3 α 3 α 3 α 4 α 4 α 4 α 5,A α 5,B α 5,C       (5.7)
The weight coefficients α i,X are defined in the range [0, 1], and their sum per phase should be equal to one. However, to allow for a small margin error the sum is constrained in the range [0.98, 1] to allow for a soft convergence of the algorithm.

     0 ≤ α i,X ≤ 1 0.98 ≤ 5 i=1 α i,X ≤ 1
(5.8)

Results

To evaluate this first approach while controlling the entire set of parameters, we use EMTP to model a complete industrial grid that represents an actual industrial grid. This model includes a section of the distribution grid, where a line fault is simulated in a parallel feeder, causing a voltage sag with a magnitude of 78% (remaining voltage) and a duration of 300 ms. The voltage and current waveforms are registered at the main monitoring point, and the three-phase voltages are used as input for the load model. The main goal is to estimate the load composition of this network by exploiting the fifth harmonic current, the active and reactive power registered during the voltage sag.

The load model is also modeled using EMTP. However, this model is a simplified version of the industrial grid, as it does not include a distribution network because the parameters and characteristics of the upstream grid are supposed to be unknown. This model is used to obtain the response to a given voltage sag profile for each load measured at the monitoring point. The injected voltages contain harmonic distortion due to the presence of power electronic loads. However, no additional upstream harmonic voltage source is considered at the busbar to model the disturbances affecting the upstream network. Both models have the same five load categories, which are described in Table 5.1. We define thirteen case scenarios with different configurations of connected loads, as presented in the table in Fig. 5.5a. The estimation errors on the weight coefficients for each scenario and for each load type are presented in the table in Fig. 5.5b. Estimation errors higher than 5% are highlighted in orange and those higher than 10% in red.

The results for eight out of thirteen cases can be considered satisfactory, with estimation errors of less than 5%. However, they are particularly high for five cases, for which we found that the sum of the individual harmonic current waveforms measured in the load model was not equal to the total harmonic current waveform measured at the monitoring point of the detailed model of the actual industrial grid. This error is not caused by the least-squares algorithm, but by the approximation made in the system of equations in (5.5). The system assumes that the load models have little impact on the voltage at the measurement point, and vise-versa. While it is the case for grids with a high level of short-circuit power, this assumption is not verified in most real cases. Although the voltage at the supply point has a limited impact on the active (P) and reactive (Q) operation point, the impact on the harmonic components can be significant because of the harmonic interactions between loads [START_REF] El-Saadany | Effect of interactions between voltage and current harmonics on the net harmonic current produced by single phase non-linear loads[END_REF].

Although partial harmonic cancellation between loads due to harmonic phase angle difference is considered when using the current waveform, the presence of multiple power electronic devices has an influence on the harmonic profile in magnitude and phase angle of the individual loads. Therefore, the responses provided by the load model are different from those obtained in the model of the complete industrial grid.

To the best of our knowledge and to this date, there are no other proposals for self-disconnected load composition estimation considering industrial loads to which we can compare our results.

Discussion, challenges, and perspectives

There are several challenges to address in order to develop an efficient algorithm for the estimation of load composition during a voltage sag, with and without load self-disconnection:

• Harmonic interaction between loads. As discussed in the previous section, harmonic interactions cause significant errors in the estimation of load composition if the total fifth harmonic current is assumed to be a simple superposition of the individual harmonic responses of the loads. Indeed, the presence of multiple non-linear loads modifies the voltage harmonic content and thus their harmonic profile, compared to when they are measured separately. However, the harmonic content of the power electronic devices provides essential information to differentiate these load categories from the others. Therefore, other estimation methods different than least-squares should be explored to better handle the non-linear coupling of harmonics.

• Influence of the supply equivalent impedance. The obtained results showed the limits of this first approach. Although the load model allows the injection of a voltage sag profile including harmonic distortion due to the presence of power electronic loads in the industrial grid, it fails to reproduce the harmonic behavior of the loads when measured individually. This could be explained by the absence of an equivalent impedance corresponding to the upstream distribution network, whose introduction in the load modeling could improve the estimation results. Unfortunately, the current version of EMTP does not allow the injection of a specific voltage profile and the inclusion of an equivalent system impedance at the same time.

• Diversity of power electronics architecture. One of the hypotheses used in this first approach was the capacity to accurately model the industrial equipment. However, for several devices, it is a real challenge to find information on their actual parameters, such as the electrical parameters of induction motors or the internal parameters of the power electronic devices, which are rarely provided by the manufacturers. Moreover, a sensitivity analysis for each load category should be performed as well, since the variation of some parameters may significantly affect the responses, particularly in terms of harmonic distortion. Estimation methods such as stochastic optimization techniques for model parameters under uncertainties could be considered. The load models could also be improved with the introduction of filters or compensation modules.

• Load aggregation. Modeling each device within the load categories in an industrial site would be tedious and complex to achieve in reality. Therefore, load aggregation techniques should be explored to reduce the total number of electrical devices into a single equivalent load model representing each load category. Some techniques have been proposed in the literature, particularly for the aggregation of induction motors [START_REF] Franklin | Improving dynamic aggregation of induction motor models[END_REF][START_REF] Kataoka | A method for aggregation of a group of induction motor loads[END_REF][START_REF] Liang | Bus-Split Algorithm for Aggregation of Induction Motors and Synchronous Motors in Dynamic Load Modeling[END_REF][START_REF] Poudel | Aggregate Model of Single Phase Induction Motors[END_REF].

• Accurate voltage tolerance curves. We have also assumed in this first approach that the voltage tolerance curves (VTC) have a perfectly rectangular shape. In reality, their shape are different, and vary according to the device's characteristics, even within the same load category [START_REF] Djokic | Sensitivity of AC Coil Contactors to Voltage Sags, Short Interruptions, and Undervoltage Transients[END_REF][START_REF] Djokic | Sensitivity of personal computers to voltage sags and short interruptions[END_REF][START_REF] Djokic | Sensitivity of AC adjustable speed drives to voltage sags and short interruptions[END_REF]. An improvement would be to integrate this diversity of VTCs and take into account their uncertainty for sensitive equipment with probability density functions, as it was proposed by Milanovic et. al in [START_REF] Milanovic | Probabilistic Assessment of Financial Losses due to Interruptions and Voltage Sags-Part I: The Methodology[END_REF].

Conclusion

Voltage sags can cause load tripping, depending on the severity of the sag and the sensitivity of industrial equipment involved in a process. It is of great interest to provide as much information as possible about the impact of a sag on the affected loads. This issue motivated the study of a methodology to estimate the self-disconnected load composition after the occurrence of a voltage sag.

Load composition estimation is one application of a broader domain called load estimation. Thus, we introduced a brief literature review on some of the most common methods for load estimation, including load modeling and parameter estimation. However, for load composition estimation, it is generally advised to use dynamic and physics-based models as they provide a more accurate response. Models based on analytical equations are not adapted for this task as the loads are aggregated into a single model regardless of their category. Similarly, static models are not well-adapted as it is the dynamic response during voltage sags that contains the most relevant information for differentiating the load categories. In general, there are fewer methods dedicated to the estimation of load composition than for stability analysis or consumption forecasting. Moreover, to the best of our knowledge, there are no other proposals in the literature that study the selfdisconnected load composition estimation after a voltage sag. This represents an opportunity to provide new solutions in the domain of voltage sag impact analysis.

Inspired by the load estimation techniques found in the literature, we proposed an approach in which we first estimate the load composition of an industrial site during a voltage sag without impact, before investigating the case of voltage sags with load self-disconnection. The results, even promising, showed some limitations of the method, particularly in terms of harmonic interactions between loads. Finally, we made several proposals to improve the load composition estimation method.

General Conclusion

The growing concern with power quality disturbances has increased these recent years due to the introduction of more sensitive and polluting power electronics devices at different levels of the electrical system. For our research, we focused on voltage sags, as they are the most frequent and impactful disturbances in industrial power grids. The financial losses linked to voltage sags are particularly high. Understanding their causes and estimating their impact are two key steps toward the implementation of adapted and cost-effective mitigation solutions.

Voltage sags are characterized in terms of amplitude, duration, phase-angle jump, and point-on-wave. The main events causing voltage sags are: line faults (balanced or unbalanced), transformer energizing, and direct motor startup. These events can take place either upstream or downstream of the monitoring point.

The core of this thesis work is dedicated to the classification of voltage sag causes and their relative location to the monitoring point. Although there is a large variety of solutions achieving high accuracy (up to 99%) for the classification of voltage sag causes, these methods encounter some limitations that prevent their application in a real industrial context, such as the requirement of large amounts of training data, low interpretability of the decision-making process and poor or unevaluated generalization capabilities across different data sources. Our objective was to propose a solution that achieves such accuracy levels while addressing these limitations.

We developed a classification algorithm based on multivariate time series signatures. The methodology follows a four-stage scheme: data acquisition, preprocessing, feature extraction, and feature analysis. It uses voltage and current waveforms as input to identify the causes of voltage sags in LV industrial grids among seven classes. The solution is based on four-dimension time series signatures, obtained through the application of the Short-Time Fourier Transform (STFT) and the Fortescue Transform. The identification stage is achieved through a distance-based classification approach referred as Nearest Neighborhood classification. Thus, the identification of a new and unknown voltage sag is achieved by comparing it to a reference signature database previously constituted. For this, we use a custom distance measure based on the dependent Dynamic Time Warping algorithm (DT W D ).

In addition, two methods for the distance-to-class calculation were proposed and compared. The most effective is based on soft-Dynamic Time Warping (soft-DTW), that is used to reduce the reference signature database into representative centroids, and in consequence significantly decreases the total computation time.

The performance of the method was analyzed in terms of class separability, prediction efficiency (accuracy and robustness to noise), and sensitivity to fundamental frequency variations. The results proved that the proposal is resilient regarding noise levels up to a SNR = 15 dB and fundamental frequency variations up to a shifting value ϵ = ±0.5Hz. Moreover, two confidence indexes (NB-KDE and RD index) were proposed and compared. The RD-index proved to be significantly more robust and stable. The information provided by this index increases the reliability of the classification process by alerting when predictions with low confidence scores are obtained and the intervention of a human expert is required, while maintaining a high degree of automation in the analysis.

Finally, the generalization capabilities of the algorithm with different data sources were evaluated using different combinations of "training" and testing data: synthetic-synthetic, real-real, synthetic-real and real-synthetic. The results were satisfactory in the four cases. The results when using a purely synthetic reference database were particularly interesting. The algorithm achieved a F1-score of 100% for the seven classes when tested on synthetic data, and a F1-score higher than 99% for the five classes existing in the database when tested with field measurement data from three different industrial sites.

In summary, the main advantages of the proposed algorithm are:

1. The reduced amount of data necessary to build the signature database (20 samples per class). 2. The possibility to constitute the signature reference database using 100% synthetic data. 3. Good generalization capabilities when implemented with field data, even if this is collected from different industrial sites. 4. Electrical interpretability of the signatures and the decision-making process. 5. Provision of a confidence index associated with the prediction.

The first three characteristics make the algorithm easy to implement in real industrial applications with no previous recorded data. The system can be developed in factory using 100% synthetic data and be directly deployed in actual industrial sites without additional training. The last two characteristics make the troubleshooting process easier and increase the general interpretability of the decision-making process, which is a requirement from industrial customers for reliability issues.

The fifth chapter focused on the impact analysis of voltage sags and, more specifically, on a methodology to estimate the self-disconnected load composition following a voltage sag. We presented a brief literature review of some of the most common methods for load estimation, including load modeling and parameter estimation methods. Most of the load estimation methods are developed for voltage stability analysis and load consumption estimation studies. The models used for these studies aim to reproduce the electrical response of aggregated loads, regardless of their individual responses. Other approaches using individual load models do not consider the devices' dynamic behavior, which is crucial for our analysis. Thus, the large majority of the methods in the literature were not adapted for the estimation of load composition.

Therefore, we proposed a first and simplified approach for this analysis. We first estimate the load composition of an industrial site for a voltage sag without impact before investigating the case of voltage sags with self-load disconnection. The results, even promising, showed some limitations, particularly in terms of harmonic interaction among the loads. Finally, we discussed some of the limits of this first approach and made several proposals to improve the load composition estimation providing guidelines for future work.

Perspectives

• Concerning the classification of voltage sag causes

The algorithm should be implemented on data collected from new industrial sites for longer periods of time to validate its effectiveness with field data on classes B1 (downstream balanced faults) and C2 (downstream transformer energizing), which where not available in the field dataset used for testing. A larger dataset could also be helpful to validate the method's performance with class C1 (upstream transformer energizing), which had less than 10 events for testing.

The methodology could also be improved to avoid classification errors in the case of voltage sags of very short duration (inferior to 30 ms), which have a risk of a peak overlap in the third ISC of their signature. Currently, these errors can be identified thanks to the confidence index. However, improvements in the feature extraction or the spatio-temporal alignment steps could be investigated.

The voltage sag cause classification algorithm is a ready-to-implement solution for low voltage (LV) sites. Nonetheless, the approach is scalable and could be extended to include new classes such as voltage sags caused by multiple faults. For instance, the algorithm could be adapted for medium voltage (MV) sites (industrial or tertiary). In MV grids, faults can cause voltage swells and a specific type of voltage sag due to single line-to-ground faults on compensated grounding systems. Preliminary observations let us think that the distinctive characteristics of these signatures should easily allow their correct identification.

Moreover, the time series classification approach presented in this PhD dissertation offers new ways to analyze short-term power quality disturbances compared to proposals in the literature. Using the entire time series reduces the risk of information loss due to the extraction of scalar features. If adequate and relevant transformations are applied, only the most meaningful characteristics of the electrical waveforms are extracted. Physically interpretable features should be preferred, since understanding their electrical behavior provides certain guarantees and reduces the amount of data necessary for the algorithm's development. In consequence, the classification stage can be significantly simplified. The transformations or number of harmonics considered in the feature extraction stage can be adapted to each case study, and the dimensions of the time series can be increased if necessary. The class separability analysis for the time series signatures proposed in this work may also be used to validate the new feature extraction process.

• Concerning the self-disconnected load composition estimation

To overcome some of the limitations of the load model, other machine learning techniques such as ANNs or DNNs could be explored. Indeed, these methods allow more flexibility to model the harmonic interactions between loads and to obtain an electrical response that also considers the short-circuit power and upstream system equivalent impedance.

Machine learning or stochastic optimization techniques could also be implemented for the parameter estimation stage. These methods would allow to integrate uncertainties in the electrical/mechanical parameters of the loads, due to the large diversity of equipment characteristics within the considered load categories.

Load aggregation techniques should be explored to reduce a group of electrical devices into a single equivalent load model representing each load category. There are some proposals in the literature, particularly for the aggregation of induction motors [START_REF] Franklin | Improving dynamic aggregation of induction motor models[END_REF][START_REF] Kataoka | A method for aggregation of a group of induction motor loads[END_REF][START_REF] Liang | Bus-Split Algorithm for Aggregation of Induction Motors and Synchronous Motors in Dynamic Load Modeling[END_REF][START_REF] Poudel | Aggregate Model of Single Phase Induction Motors[END_REF]. Finally, more realistic representations of the VTCs should be investigated, for instance, with probability density functions as proposed in [START_REF] Milanovic | Probabilistic Assessment of Financial Losses due to Interruptions and Voltage Sags-Part I: The Methodology[END_REF].

Résumé en français

L'augmentation de la consommation d'électricité, les nouveaux usages et la nécessité de préserver l'environnement sont à l'origine de la transition énergétique. Cette transition modifie profondément les réseaux électriques avec une plus grande pénétration des énergies renouvelables dans le mix énergétique, à la fois associée à des interfaces d'électronique de puissance et à la digitalisation du contrôle-commande [START_REF] Steimer | Enabled by high power electronics -Energy efficiency, renewables and smart grids[END_REF][START_REF] Tang | Power electronics: The enabling technology for renewable energy integration[END_REF]. Les flux d'énergie deviennent multidirectionnels et nécessitent une gestion plus intelligente pour répondre aux exigences techniques en termes de disponibilité, fiabilité, sécurité et qualité de l'alimentation électrique. L'analyse de la qualité d'électricité est ainsi devenue une préoccupation en forte hausse pour les fournisseurs d'énergie de même que pour leurs clients au cours des dernières années. Un approvisionnement énergétique fiable garantit le fonctionnement optimal des équipements électriques du réseau, alors qu'une mauvaise qualité d'électricité peut avoir comme conséquence l'interruption des lignes de production ou des services, le dysfonctionnement des équipements, et leur endommagement allant jusqu'à la casse. Les pertes financières associées peuvent être importantes [START_REF] Targosz | Pan-European power quality survey[END_REF] pour les industriels et les grands clients tertiaires tels que des hôpitaux, des datacenters, etc. Ainsi, identifier l'origine des perturbations de la qualité de l'énergie et évaluer leurs effets sur les équipements industriels est essentiel pour identifier et proposer des solutions palliatives qui soient à la fois adaptées et rentables, afin de réduire l'impact sur la productivité du site. Aujourd'hui, l'analyse des perturbations de la qualité d'électricité est le plus souvent réalisée par des experts du domaine. Elle nécessite un haut niveau de connaissances et d'expertise pour établir un diagnostic fiable, et proposer des solutions pertinentes. Cependant, ce processus est très chronophage et nécessite une intervention spécifique sur site pour l'acquisition des données. De plus, certaines perturbations électriques sont rares en termes de localisation et de fréquence. Leur enregistrement peut donc prolonger encore l'étape d'acquisition des données. Des dispositifs de surveillance tels que des qualimètres peuvent alors être placés en permanence pour faire face à ce problème. Néanmoins, le traitement et l'analyse de grandes quantités de données peuvent également prendre un temps conséquent.

De ce fait, ce travail de recherche vise à développer un système intelligent d'analyse des perturbations en traitant les mesures électriques (tension et courant triphasés) provenant d'un seul dispositif de mesure placé au point d'alimentation électrique principal d'un site industriel ou tertiaire.

Parmi les principales perturbations affectant les réseaux industriels, les creux de tension sont les plus fréquents et les plus impactants [START_REF]IEEE Recommended Practice for Monitoring Electric Power Quality[END_REF][START_REF] Wagner | Power quality and factory automation[END_REF][START_REF] Sarmiento | A voltage sag study in an industry with adjustable speed drives[END_REF], car ils peuvent provoquer le dysfonctionnement de certains équipements ainsi que l'arrêt intempestif des processus industriels. Pour ces raisons, nous nous sommes focalisés sur cette type de perturbation. L'analyse des creux de tension s'effectue en deux étapes: classification de sources de creux de tension et analyse de leur impact sur les équipements industriels.

Cependant, il y a un certain nombre de défis à rélever afin de développer une telle méthodologie. Parmi les plus importants nous pouvons noter :

• L'accès à un unique point de mesure au niveau du point de raccordement principal du site. L'objectif du système de diagnostic est d'être le moins intrusif possible, avec des données collectées à partir d'un seul point de mesure.

• Le manque d'informations sur la topologie du réseau industriel (interne) ainsi que sur les équipements industriels en aval du point de mesure.

• La diversité des appareils électriques, avec différents niveaux de sensibilité aux creux de tension. La tolérance aux creux de tension dépend non seulement des caractéristiques des creux, mais aussi des caractéristiques de l'équipement industriel. Le développement d'un système capable de s'adapter à cette diversité d'équipements répresente un vrai défi.

Dans le Chapitre I, nous présentons le contexte électrique nécessaire à la compréhension et au développement de méthodes d'analyse automatique des perturbations électriques. Une attention particulière est portée aux creux de tension. Après la caractérisation des creux de tension en termes d'amplitude, de durée, de saut de phase et de point sur l'onde, nous abordons leurs principales causes: défauts de ligne, enclenchement de transformateurs et démarrage direct de moteurs asynchrones. De la même manière, nous décrivons les principaux effets et conséquences des creux de tension sur les équipements industriels les plus sensibles: moteurs à induction, variateurs de vitesse, ordinateurs et PLC, contacteurs AC et lampes à décharge. Les principales causes et conséquences de la distorsion harmonique sont également abordés. La Figure 1 présente un schéma simplifé d'un réseau électrique industriel en basse tension (BT). Dans le Chapitre II, une étude bibliographique dans le domaine de l'analyse de la qualité d'électricité est présentée, avec un focus sur la classification de perturbations. Nous présentons ensuite un comparatif entre les différentes méthodes dans la littérature pour la classification des causes creux de tension. Bien que la précision de classification globale des méthodes soit élevée (entre 92% et 100%), nous avons identifié certaines de leurs principales limitations, afin de les prendre en compte dans le développement de notre solution, parmi lesquelles les plus importantes sont :

• La nécessité d'avoir accès à d'importantes quantités de données pour le développement des algorithmes, notamment dû à l'utilisation des méthodes statistiques pour la grande majorité des solutions dans la litérature. Ces méthodes sont fortement dépendantes de données d'entraînement. L'accès à ces données étant difficile, la quantité minimale requise doit être limitée. Notre objectif est de proposer une solution qui atteigne de tels standards, mais qui en même temps traite les limitations mentionnées ci-dessus. Nous tenons en compte également que l'algorithme à développer a comme finalité d'être deployé et commercialisé. Sa mise en oeuvre se doit donc d'être réaliste d'un point de vue industriel.

Dans le Chapitre III, nous présentons une nouvelle méthodologie de classification des causes des creux de tension. Les sept classes définies sont les suivantes : défaut amont équilibré, défaut amont déséquilibré, défaut aval équlibré, défaut aval déséquilibré, enclenchement transformateur amont, enclenchement transformateur aval et démarage moteur aval. La solution proposée est basée sur l'extraction d'indicateurs interprétables, intégrant des connaissances expert. Ces indicateurs sont des signatures de séries temporelles multivariées qui, comparées aux indicateurs scalaires, présentent un risque plus faible de perte d'information. Ensuite, nous proposons une approche par classification des séries temporelles pour classer les signatures. Bien que cette méthode soit basée sur une base de données d'entraîenement (référence), nous démontrons par la suite que la quantité de données requise pour le développement de la solution est bien inférieure à celle des autres approches de la littérature. En outre, nous nous assurons que les données requises soient accessibles par simulation et que l'algorithme développé soit capable de fonctionner efficacement dans différents types de réseaux industriels. Comme illustré sur la Figure 2, la méthode est composée de quatre étapes: Les résultats prouvent que la méthode est résiliente face à des niveaux de bruit allant jusqu'à un niveau SNR = 15 dB et des variations de fréquence fondamentale jusqu'à ϵ = ±0, 5 Hz. L'algorithme a été évalué sur de données croisées, synthétiques et réelles, montrant de très bonnes capacités de généralisation. L'algorithme atteint un F1-score de 100% sur les sept classes définies avec une base de données synthétique réduite (140 événements) lorsqu'il a été testé sur des données synthétiques. Cependant, les résultats les plus intéressans ont été obtenus avec une base de données de signatures synthétiques lors de l'entraînement et appliqué à des données de terrain provenant de trois sites industriels différents pour la partie test. L'algorithme atteint ici un F1-score de 99,32% pour cinq des sept classes définies.

Le Chapitre V traite sur l'analyse de l'impact des creux de tension et, plus particulièrement, sur une méthodologie permettant d'estimer la composition de la charge auto-déconnectée suite au creux. La déconnexion de charges étant en fonction de la gravité du creux et de la sensibilité des équipements industriels.

Il est important de souligner qu'à notre connaissance, il n'existe pas à ce jour de propositions dans la littérature qui étudient l'estimation de la composition de la charge auto-déconnectée après un creux de tension. Pour cette raison, nous présentons une brève étude bibliographique sur un sujet plus large qui est celui de l'estimation de la charge ou load estimation en anglais.

Ainsi, un algorithme de classification basé sur des signatures de séries temporelles multivariées a été développé. La méthodologie suit un schéma en quatre étapes afin d'identifier les causes de creux de tension dans les réseaux industriels BT parmi les sept classes définies. Les principaux avantages de cette méthodologie sont :

• La quantité réduite de données nécessaires pour construire la base de données de signatures de référence, qui peut être entièrement composée de données synthétiques.

• La fourniture d'un indice de confiance associé à la prédiction.

• L'interprétabilité électrique des signatures et le processus de décision.

• La robustesse aux niveaux de bruit jusqu'à un SNR = 15 dB et aux variations de la fréquence fondamentale jusqu'à ±0.5Hz.

• Ses bonnes capacités de généralisation lorsque l'algorithme est appliqué sur de données de terrain réelles, même pour des sites industriels différents.

Les trois premières caractéristiques rendent l'algorithme facile à mettre en oeuvre dans des applications industrielles réelles. Le système peut être développé en amont en utilisant 100% de données synthétiques et être directement déployé sur des sites industriels réels. Les deux dernières caractéristiques facilitent le processus d'analyse d'erreurs et augmentent l'interprétabilité générale de la solution.

La deuxième partie de la thèse a été consacrée à l'analyse de l'impact des creux de tension et, plus précisément, à une méthodologie d'estimation de la composition de la charge auto-déconnectée suite à un creux de tension. Une première approche pour cette analyse a été proposée. Nous avons d'abord estimé la composition de la charge d'un site industriel pour une chute de tension sans impact avant d'étudier le cas des chutes de tension avec déconnexion de charges. Les résultats, même prometteurs, ont montré certaines limites, notamment lorsqu'il existe de l'interaction harmonique entre équipements. Enfin, nous avons analysé les limites de cette première approche et nous avons fait plusieurs propositions pour améliorer la méthode en fournissant des pistes pour de travaux futurs.

En ce qui concerne les perspectives elles peuvent être envisagées dans les deux directions suivantes :

• Concernant la classification des causes des creux de tension 1. L'algorithme peut être appliqué sur des données provenant de nouveaux sites industriels pendant des périodes plus longues afin de valider son efficacité sur les deux classes manquantes.

2. La méthodologie pourrait également être améliorée pour éviter les erreurs de classification dans le cas de creux de tension de très courte durée (inférieures à 30 ms), qui présentent un risque d'erreur de classification.

3. L'approche étant évolutive, elle pourrait être étendue pour inclure de nouvelles classes telles que les creux de tension causés par des défauts évolutifs.

De plus, l'approche par classification des séries temporelles présentée dans ces travaux offre de nouvelles façons d'analyser les perturbations électriques courtterme par rapport aux méthodes dans la littérature. L'utilisation intégrale des séries temporelles permettant de réduire le risque de perte d'information due à l'extraction d'indicateurs scalaires.

1 . 1 107 xi 4 . 3 4 . 4 4 . 6 4 . 7 4 . 9 137 5. 3 143 xiiiof Tables 1 . 1 4 . 3 PFC 1 V 1 V 2 VFNQ

 11107434446474913731431143112 Simplified diagram of a distribution network. The monitoring device at a LV industrial site is illustrated. . . . . . . . . . . . . . . . . . . 1.2 Different voltage sag representations . . . . . . . . . . . . . . . . . 1.3 Voltage sag representations of a double-phase fault . . . . . . . . . 1.4 Voltage sag representations due to a transformer energizing . . . . . 1.5 Voltage sag representations due to an induction motor direct startup 1.6 CBEMA, ITIC and SEMI F47 curves, from [1] . . . . . . . . . . . 1.7 Voltage-tolerance curve (VTC) . . . . . . . . . . . . . . . . . . . . . 1.8 Harmonic distortion of current and decomposition in harmonic components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Power quality analysis main categories . . . . . . . . . . . . . . . . 2.2 Generalized flowchart for power quality disturbance classification . . 2.3 RMS values of a synthetic voltage sag generated with a numerical model (a) and a simulation model (b). Only the simulation model emulates the influence of the induction motor (IM) on the sag. . . . 2.4 Format transformation from voltage waveforms to (a) RMS values and (b) and gray-scale images, for voltage sag source identification based on image processing techniques. [2]. . . . . . . . . . . . . . . 3.1 Flowchart of the two modes: reference database constitution and classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 (a) Centralized and (b) decentralized approaches for data collection and processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.3 Simplified diagram of the simulated industrial network. . . . . . . . ix LIST OF FIGURES 3.4 Short-Time Fourier transform calculation . . . . . . . . . . . . . . . 3.5 Influence of window type on the STFT result . . . . . . . . . . . . . 3.6 Influence of window hop length W H on the STFT result . . . . . . . 3.7 One-phase current signal (a) and the module of its corresponding Short-Time Fourier Transform (b). . . . . . . . . . . . . . . . . . . 3.8 Fortescue transform applied to a three-phase unbalanced system. . . 3.9 Fortescue transform applied to three-phase voltage time series (harmonic 1 or fundamental), resulting in three instantaneous symmetrical components . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Waveform downsampling effect on the extracted signatures . . . . . 3.11 RMSE (%) as a function of the sampling frequency of the electrical waveforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Feature extraction process and the obtained components. . . . . . . 3.13 Voltage sag caused by an upstream unbalanced fault (field data). . . 3.14 Voltage sag caused by an upstream transformer energizing (field data). 3.15 Voltage sag caused by a downstream induction motor startup (field data). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 Two voltage sags caused by upstream unbalanced faults (A2) . . . . 3.17Voltage positive-sequence harmonic 1: (a) fault, (b) transformer energizing and (c) motor startup . . . . . . . . . . . . . . . . . . . 3.18 Voltage negative-sequence harmonic 1: (a) unbalanced fault, (b) balanced fault, (c) transformer energizing and (d) motor startup . . 3.19 Voltage positive-sequence harmonic 2: (a) fault, (b) transformer energizing, (c) motor direct startup . . . . . . . . . . . . . . . . . . 3.20 Current positive-sequence harmonic 1: (a) downstream fault, (b) downstream transformer energizing, (c) downstream motor startup (d) upstream fault and (e) upstream transformer energizing . . . . 3.21 Time series matching using (a) Euclidean distance and (b) Dynamic Time Warping algorithm. . . . . . . . . . . . . . . . . . . . . . . . . 3.22 Time alignment of two univariate time series. The optimal warping path is given by the lowest value of the cumulative cost E, obtained from the distance matrix A. . . . . . . . . . . . . . . . . . . . . . . x 3.23 Time alignment of (a) two multivariate time series with (b) Dependent DTW and (c) Independent DTW [3]. . . . . . . . . . . . . . . 85 3.24 Poor time alignment due to singularity points (red) can lead to high distortion of time series. . . . . . . . . . . . . . . . . . . . . . . . . 86 3.25 Global constraints on the warping path (a) Sakoe-Chiba band and (b) Itakura parallelogram . . . . . . . . . . . . . . . . . . . . . . . . 87 3.26 Local step patterns: (a) symmetric2, (b) symmetric1, (c) asymmetric 87 3.27 Time alignment results of two signatures of the same class, using (a) symmetric2 and (b) symmetric1 step patterns. . . . . . . . . . . 89 3.28 Time alignment results of two signatures belonging to different classes, using (a) symmetric2 and (b) symmetric1 step patterns. . . . . . . 90 3.29 Time alignment results of two signatures belonging to the same class, using (a) asymmetric and (b) symmetric1 step patterns. . . . 91 3.30 Result of the space and time alignment step. . . . . . . . . . . . . . 92 3.31 1-Nearest Neighbor classifier using a custom distance. . . . . . . . . 93 3.32 Calculation of the distance to a class C k using bootstrapping . . . . 95 3.33 Nearest Neighborhood with centroid estimation using soft-DTW. . . 96 3.34 Centroid estimation for multivariate time series using (a) a standard Euclidean averaging, (b) DBA with expectation minimization, (c) soft-Dynamic Time Warping (γ = 0.5). . . . . . . . . . . . . . . . . 98 3.35 Centroid estimation using soft-DTW with different values of γ. . . . 99 3.36 Optimization of the parameter γ . . . . . . . . . . . . . . . . . . . . 100 3.37 (a) Likelihood function and (b) output score of a binary Naive Bayes classifier using KDE for a new voltage sag signature q * to the class C 2 , given the distance to the class x = d 2 (q * ). . . . . . . . . . . . . 102 3.38 (a) Distances d k (q * ) of a new sag q * to each class k, and the corresponding (b) membership probabilities p k (q * ). . . . . . . . . . . . . 103 4.1 Intra-class and inter-class distances . . . . . . . . . . . . . . . . . . 106 4.2 Intra-class distribution H C2 and inter-class distribution H C2-C1 , with BC = 0.07. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bootstrapping method: (a) F1-score and mean values for NB-KDE and RD indexes and the corresponding (b) Log-loss values associated to the output predictions for 16 different noise levels. . . . . . Centroid estimation method: (a) F1-score and mean values for NB-KDE and RD indexes and the corresponding (b) Log-loss values associated to the output predictions for 16 different noise levels. . . 4.5 Varying noise levels on the voltage waveform and their impact to the signatures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bootstrapping method: F1-score and mean values for NB-KDE and RD indexes (a) and the corresponding Log-loss values (b) associated to the output predictions, for 5 different levels of fundamental frequency variation. . . . . . . . . . . . . . . . . . . . . . . . . . . . Centroid estimation method: F1-score and mean values for NB-KDE and RD indexes (a) and the corresponding Log-loss values (b) associated to the output predictions, for 5 different levels of fundamental frequency variation. . . . . . . . . . . . . . . . . . . . 4.8 Effects of the fundamental frequency variations on the second ISC of a voltage sag caused by a transformer energizing. . . . . . . . . . Bootstrapping method: (a) F1-score and mean values for NB-KDE and RD indexes and the corresponding (b) Log-Loss values associated to the output predictions, for different sizes of the synthetic database. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Centroid estimation method: (a) F1-score and mean values for NB-KDE and RD indexes and the corresponding (b) Log-Loss values associated to the output predictions, for different sizes of the synthetic database. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Classification errors in experiment no. 3. Reference signature of class A1 (50 ms) and query signature of class A2 (30 ms) . . . . . . 4.12 Classification errors in experiment no. 3. Reference signature of class A2 (50 ms) and query signature of class A2 (30 ms) . . . . . . 4.13 Analysis of misclassified events for experiment 3: (a) Boxplot of the RD confidence index, 5 events raise an alert due to low confidence index, with two of them corresponding to the misclassification errors and (b) detailed membership probabilities according to RD index for the events raising an alert. . . . . . . . . . . . . . . . . . . . . . xii 5.1 (a) Voltage sag without impact and (b) voltage sag with load selfdisconnection or load tripping . . . . . . . . . . . . . . . . . . . . . 129 5.2 Two-step approach for load composition estimation after a voltage sag with and without impact . . . . . . . . . . . . . . . . . . . . . . Diagram of the load model used for the load composition estimation during a voltage sag . . . . . . . . . . . . . . . . . . . . . . . . . . 138 5.4 Measurement-based approach for load composition estimation . . . 140 5.5 (a) Description of cases with different configurations of connected loads and (b) estimation errors per case and per load . . . . . . . . List Typical characteristics of Power Quality Disturbances . . . . . . . . 1.2 Mixed networks: voltage sag incidence -95th percentile [4] . . . . . 1.3 Current waveforms and harmonic content of non-linear devices . . . 1.4 Effects of current harmonics on the calculation of power quantities. Voltage harmonics are considered negligible. . . . . . . . . . . . . . 2.1 Numerical modeling of power quality disturbances . . . . . . . . . . 2.2 Benefits and drawbacks of signal processing techniques . . . . . . . 2.3 Example of common extracted features using different signal processing techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Main properties of classifiers used for power quality classification . . 2.5 Comparative review of Voltage Sag Cause Classification methods in the literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Definition of classes . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Real dataset description . . . . . . . . . . . . . . . . . . . . . . . . 3.3 LV industrial loads . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Step pattern comparison . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Comparison between centroid estimation methods . . . . . . . . . . 4.1 Bhattacharyya coefficient (BC) between the intra-class and interclass distance distributions . . . . . . . . . . . . . . . . . . . . . . . 4.2 Prediction and confidence index results using synthetic data -Mean distance with Bootstrapping method . . . . . . . . . . . . . . . . . xiv Prediction and confidence index results using synthetic data -Centroid estimation method . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Computational cost comparison . . . . . . . . . . . . . . . . . . . . 4.5 Description of the experiments . . . . . . . . . . . . . . . . . . . . . 4.6 Experiment 1: Training and testing with synthetic data . . . . . . . 4.7 Experiment 2: Training and testing with real field data . . . . . . . 4.8 Experiment 3: Training with synthetic data, and testing with real field data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Experiment 4: Training with real field data, and testing with synthetic data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Comparison between two distance-to-class calculation methods and two confidence scores . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Load composition estimation -Load description . . . . . . . . . . . xv Glossary List of Abbreviations Dy Delta-wye winding connection 1NN-DTW 1-Nearest Neighbors with Dynamic Time Warping AC Alternative current ANN Artificial Neural Network ASD Adjustable speed drive CBEMA Computer and Business Equipment Manufacturer'Naives Bayes using a Kernel Density Estimator, confidence index NMSE Normalized Mean Squared Error NRMSE Normalized Root Mean Squared Error PCC Point of common coupling r) Cost matrix (DTW) γ Smoothing parameter (soft-DTW) γ opt Optmimal value for smoothing parameter (soft-DTW) rk Barycenter or centroid of class C k A Distance matrix D(q, r) Distance between two time series signatures q and r H α-β Inter-class distance distribution of class α and class β H α Intra-class distance distribution of class α ϕ(l) Warping path I + {h1} Current positive-sequence harmonic + {h1} Voltage positive-sequence harmonic + {h2} Voltage positive-sequence harmonic -{h1} Voltage negative-sequence harmonic 1 BC Bhattacharyya coefficient C k Any class from the defined list of classes D Dimension of the time series signature D Distortion power (KVARD) d k (q * ) Distance between signature q * to class C k Nominal Total number of reference signatures in the database (all classes) N k Size of the class C k P Active power (kW) p k (q * ) Membership probability of signature q * to the class C k Reactive Reference signatures from class C k S Reactive power (KVA) Z (%) Distortion rate xix

Figure 1 . 1 :

 11 Figure 1.1: Simplified diagram of a distribution network. The monitoring device at a LV industrial site is illustrated.

Figure 1 . 2 :

 12 Figure 1.2: Different voltage sag representations

Figure 1 . 3 :

 13 Figure 1.3: Voltage sag representations of a double-phase fault

Fig. 1 .

 1 Fig. 1.4 illustrates a voltage sag at 0.05s caused by a transformer energizing. The exponential shape of the voltage recovery and the voltage unbalance are clearly visible in the RMS plot.

Figure 1 . 4 :

 14 Figure 1.4: Voltage sag representations due to a transformer energizing

Figure 1 . 5 :

 15 Figure 1.5: Voltage sag representations due to an induction motor direct startup

Figure 1 . 6 :

 16 Figure 1.6: CBEMA, ITIC and SEMI F47 curves, from [1]

Figure 1 . 8 :

 18 Figure 1.8: Harmonic distortion of current and decomposition in harmonic components

Figure 2 . 1 :

 21 Figure 2.1: Power quality analysis main categories

Figure 2 . 2 :

 22 Figure 2.2: Generalized flowchart for power quality disturbance classification

Figure 2 . 3 :

 23 Figure 2.3: RMS values of a synthetic voltage sag generated with a numerical model (a) and a simulation model (b). Only the simulation model emulates the influence of the induction motor (IM) on the sag.

  (a) RMS values of voltage sags (b) Gray-scale images of voltage sags

Figure 2 . 4 :

 24 Figure 2.4: Format transformation from voltage waveforms to (a) RMS values and (b) and gray-scale images, for voltage sag source identification based on image processing techniques. [2].

Figure 3 . 2 :

 32 Figure 3.2: (a) Centralized and (b) decentralized approaches for data collection and processing

Figure 3 . 3 :

 33 Figure 3.3: Simplified diagram of the simulated industrial network.

  Fig. 3.4 illustrates the calculation process of the STFT.

Figure 3 . 4 :

 34 Figure 3.4: Short-Time Fourier transform calculation

Figure 3 . 5 :Figure 3 . 6 :Figure 3 . 7 :

 353637 Figure 3.5: Influence of window type on the STFT result

  .8 and described in equation (3.3).

Figure 3 . 8 :

 38 Figure 3.8: Fortescue transform applied to a three-phase unbalanced system.

Figure 3 . 9 :

 39 Figure 3.9: Fortescue transform applied to three-phase voltage time series (harmonic 1 or fundamental), resulting in three instantaneous symmetrical components

Figure 3 . 10 :Figure 3 . 11 :

 310311 Figure 3.10: Waveform downsampling effect on the extracted signatures

Figure 3 . 12 :

 312 Figure 3.12: Feature extraction process and the obtained components.

Figure 3 . 13 :

 313 Figure 3.13: Voltage sag caused by an upstream unbalanced fault (field data).

Figure 3 . 14 :

 314 Figure 3.14: Voltage sag caused by an upstream transformer energizing (field data).

Figure 3 . 15 :

 315 Figure 3.15: Voltage sag caused by a downstream induction motor startup (field data).

Figure 3 . 16 :

 316 Figure 3.16: Two voltage sags caused by upstream unbalanced faults (A2)

Figure 3 . 17 :

 317 Figure 3.17: Voltage positive-sequence harmonic 1: (a) fault, (b) transformer energizing and (c) motor startup

Figure 3 . 18 : 3 .

 3183 Figure 3.18: Voltage negative-sequence harmonic 1: (a) unbalanced fault, (b) balanced fault, (c) transformer energizing and (d) motor startup

Figure 3 . 19 : 4 .

 3194 Figure 3.19: Voltage positive-sequence harmonic 2: (a) fault, (b) transformer energizing, (c) motor direct startup

Figure 3 . 20 :

 320 Figure 3.20: Current positive-sequence harmonic 1: (a) downstream fault, (b) downstream transformer energizing, (c) downstream motor startup (d) upstream fault and (e) upstream transformer energizing

Figure 3 . 21 :

 321 Figure 3.21: Time series matching using (a) Euclidean distance and (b) Dynamic Time Warping algorithm.

. 8 )Figure 3 . 22 :

 8322 Figure 3.22: Time alignment of two univariate time series. The optimal warping path is given by the lowest value of the cumulative cost E, obtained from the distance matrix A.

Figure 3 . 23 :

 323 Figure 3.23: Time alignment of (a) two multivariate time series with (b) Dependent DTW and (c) Independent DTW [3].

Figure 3 . 24 :

 324 Figure 3.24: Poor time alignment due to singularity points (red) can lead to high distortion of time series.

Figure 3 . 25 :

 325 Figure 3.25: Global constraints on the warping path (a) Sakoe-Chiba band and (b) Itakura parallelogram

Figure 3 . 26 :

 326 Figure 3.26: Local step patterns: (a) symmetric2, (b) symmetric1, (c) asymmetric

( a )

 a Symmetric2 step pattern (b) Symmetric1 step pattern

Figure 3 . 27 :

 327 Figure 3.27: Time alignment results of two signatures of the same class, using (a) symmetric2 and (b) symmetric1 step patterns.

( a )

 a Symmetric2 step pattern (b) Symmetric1 step pattern

Figure 3 . 28 :

 328 Figure 3.28: Time alignment results of two signatures belonging to different classes, using (a) symmetric2 and (b) symmetric1 step patterns.

( a )

 a Asymmetric step pattern (b) Symmetric1 step pattern

Figure 3 . 29 :

 329 Figure 3.29: Time alignment results of two signatures belonging to the same class, using (a) asymmetric and (b) symmetric1 step patterns.

Fig. 3 .

 3 Fig. 3.30 illustrates the result of the space and time-alignment step using the symmetric1 step pattern between two voltage sags caused by upstream transformer energizing.

  (a) Before space and time-alignment (b) After space and time-alignment

Figure 3 . 30 :

 330 Figure 3.30: Result of the space and time alignment step.

Figure 3 .

 3 Figure 3.31: 1-Nearest Neighbor classifier using a custom distance.

Figure 3 . 32 :

 332 Figure 3.32: Calculation of the distance to a class C k using bootstrapping

Figure 3 . 33 :

 333 Figure 3.33: Nearest Neighborhood with centroid estimation using soft-DTW.

Figure 3 . 34 :

 334 Figure 3.34: Centroid estimation for multivariate time series using (a) a standard Euclidean averaging, (b) DBA with expectation minimization, (c) soft-Dynamic Time Warping (γ = 0.5).

  (a) γ = 0.01 (b) γ = 10

Figure 3 . 35 :

 335 Figure 3.35: Centroid estimation using soft-DTW with different values of γ.

Figure 3 . 36 :

 336 Figure 3.36: Optimization of the parameter γ

Figure 3 .

 3 Figure 3.37: (a) Likelihood function and (b) output score of a binary Naive Bayes classifier using KDE for a new voltage sag signature q * to the class C 2 , given the distance to the class x = d 2 (q * ).

Fig. 3 .

 3 Fig. 3.38 illustrates the distance-to-class and the corresponding membership probabilities of a new voltage sag caused by an upstream balanced fault (class

  (a) Distance-to-class (b) Membership probability

Figure 3 .

 3 Figure 3.38: (a) Distances d k (q * ) of a new sag q * to each class k, and the corresponding (b) membership probabilities p k (q * ).

Fig. 4 .

 4 1 illustrates both intra-class and inter-class distances.

Figure 4 . 1 :

 41 Figure 4.1: Intra-class and inter-class distances

Fig. 4 .

 4 2 illustrates the intra-class distribution of class C2 (downstream transformer energizing) and the inter-class distribution between C2 and C1 (upstream transformer energizing) 1 .

Figure 4 . 2 :

 42 Figure 4.2: Intra-class distribution H C2 and inter-class distribution H C2-C1 , with BC = 0.07.

Figure 4 . 3 :

 43 Figure 4.3: Bootstrapping method: (a) F1-score and mean values for NB-KDE and RD indexes and the corresponding (b) Log-loss values associated to the output predictions for 16 different noise levels.

Figure 4 . 4 :

 44 Figure 4.4: Centroid estimation method: (a) F1-score and mean values for NB-KDE and RD indexes and the corresponding (b) Log-loss values associated to the output predictions for 16 different noise levels.

Figure 4 . 5 :

 45 Figure 4.5: Varying noise levels on the voltage waveform and their impact to the signatures

  (a) F1-score and mean values for NB-KDE and RD confidence indexes (b) Log-loss values for NB-KDE and RD confidence indexes

Figure 4 . 6 :

 46 Figure 4.6: Bootstrapping method: F1-score and mean values for NB-KDE and RD indexes (a) and the corresponding Log-loss values (b) associated to the output predictions, for 5 different levels of fundamental frequency variation.

  (a) F1-score and mean values for NB-KDE and RD confidence indexes (b) Log-loss values for NB-KDE and RD confidence indexes

Figure 4 . 7 :

 47 Figure 4.7: Centroid estimation method: F1-score and mean values for NB-KDE and RD indexes (a) and the corresponding Log-loss values (b) associated to the output predictions, for 5 different levels of fundamental frequency variation.

Figure 4 . 8 :

 48 Figure 4.8: Effects of the fundamental frequency variations on the second ISC of a voltage sag caused by a transformer energizing.

Figure 4 . 9 :

 49 Figure 4.9: Bootstrapping method: (a) F1-score and mean values for NB-KDE and RD indexes and the corresponding (b) Log-Loss values associated to the output predictions, for different sizes of the synthetic database.

Figure 4 . 10 :

 410 Figure 4.10: Centroid estimation method: (a) F1-score and mean values for NB-KDE and RD indexes and the corresponding (b) Log-Loss values associated to the output predictions, for different sizes of the synthetic database.

Figure 4 . 11 :

 411 Figure 4.11: Classification errors in experiment no. 3. Reference signature of class A1 (50 ms) and query signature of class A2 (30 ms)

Figure 4 . 12 :

 412 Figure 4.12: Classification errors in experiment no. 3. Reference signature of class A2 (50 ms) and query signature of class A2 (30 ms)

Figure 4 . 13 :

 413 Figure 4.13: Analysis of misclassified events for experiment 3: (a) Boxplot of the RD confidence index, 5 events raise an alert due to low confidence index, with two of them corresponding to the misclassification errors and (b) detailed membership probabilities according to RD index for the events raising an alert.

Figure 5 . 1 :Fig. 5

 515 Figure 5.1: (a) Voltage sag without impact and (b) voltage sag with load selfdisconnection or load tripping

  Fig. 5.3 illustrates the approach.

Figure 5 . 3 :

 53 Figure 5.3: Diagram of the load model used for the load composition estimation during a voltage sag

Figure 1 :

 1 Figure 1: Schéma simplifié d'un réseau électrique industriel en BT.

Figure 2 :

 2 Figure 2: Schéma global de l'algorithme avec deux modes de fonctionnement: constitution de la base de données et classification d'un nouveau creux de tension.

Figure 3 :

 3 Figure 3: Alignement spatio-temporel entre deux signatures appartenant à la même classe.

Figure 4 :

 4 Figure 4: Analyse d'erreurs de classification en utilisant l'indice RD, où 5 événements déclenchent une alerte en raison d'un faible indice de confiance, et dont deux correspondent à des erreurs de classification.

Table 1 .

 1 1: Typical characteristics of Power Quality Disturbances

	Disturbance	Duration	Spectral content Amplitude Cause		Consequence
					Switching event,
	Transient (im-pulsive, oscilla-tory)	<10 ms	5 kHz -5 MHz	0-8 pu	transformer or capacitor bank energizing, light-	Equipment dys-function
					ning
	Voltage sag	10 ms -1 min		0.1-0.9 pu	Line motor startup, faults, transformer energizing	Motor stalling, equipment's protection trip-ping, stopped processes
							Equipment's
					Line faults, large	protection trip-
	Voltage swell	10 ms -1 min		1.1-1.2 pu	load switching,	ping,	overcur-
					load shedding	rent, equipment
							damage
	Short interrup-tion	10 ms -1 min		<0.1 pu	Line equipment mal-faults, function, loose connections	Equipment's protection trip-ping, overcur-rent, equipment damage
							Heating,	pre-
					Non-linear loads	mature	aging,
	Harmonics	steady-state	0-2 kHz 1	0-20%	(power tronic devices, elec-	transformer saturation,
					arc furnace)	performance
							reduction
					Non-linear and
					fluctuating loads
	Interharmonics	steady-state	0-2 kHz 1	0-2%	(arc welding	furnace, ma-	Light flicker
					chine,	wind
					turbine)
	Imbalance	steady-state		0.5-5%	Asymmetrical loads, one-phase loads	Motor heating, vibrations
					Arc	furnace,
	Flicker	steady-state		0.1-7%	welding chines,	ma-wind	Light flicker
					turbines
					Phase	current	Frequency	or
	Notching	steady-state			commutation in power electron-	timing errors in power electron-
					ics		ics

Table 1 .

 1 2 presents the results of the UNIPEDE DISDIP measurement campaign[START_REF]Voltage dips and short interruptions in medium voltage public electricity supply systems[END_REF] for 52 measurement sites with mixed networks (overhead and underground). Nine countries of Europe co-operated to determine the number of voltage sags and short interruptions per year experienced by electricity clients connected to LV and MV networks.

Table 1 .

 1 2: Mixed networks: voltage sag incidence -95th percentile[START_REF]Voltage dips and short interruptions in medium voltage public electricity supply systems[END_REF]

	Duration t

Table 1 .

 1

	3: Current waveforms and harmonic content of non-linear devices
	Equipment	Current waveform	Harmonic spectrum
	Single-phase
	rectifier
	Three-phase
	rectifier (6-pulse bridge)
	LED lamp
	Transformer energizing
	(after 150 ms)

Table 1 .

 1 4: Effects of current harmonics on the calculation of power quantities. Voltage harmonics are considered negligible.

Table 2 .

 2 3 presents a brief description of some examples of features extracted based on different signal processing tools.

Table 2 .

 2 2: Benefits and drawbacks of signal processing techniques

	SP technique	Benefits			Drawbacks		Ref.
	Short-time Fourier Trans-	Useful for time-frequency	Limited time-frequency	[37,	68,
	form (STFT)	analysis of non-stationary	resolution due to fixed	69]
		signals. Easy to imple-	window.
		ment.
	Wavelet Transform (WT) Allows a fine decomposi-	Performance	depends	[71,	72,
		tion of the spectral con-	on the choice of mother	73, 74, 75,
		tent in time.			wavelet.	High level of	76, 77, 78]
					decomposition can be
					computationally expen-
					sive. Very noise sensitive.
	S-Transform (ST)	Suitable	for	time-	Not suitable for harmonic	[81,	66,
		frequency		analysis.	analysis	because	the	82, 83, 84,
		Good trade-off between	widths of frequency win-	85]
		time and frequency res-	dows are directly related
		olution thanks to its	to their central frequency.
		frequency-dependent vari-	High computational cost.
		able window.
	Hilbert-Huang Transform	Able to identify subtle	Limited	performance	[87,	88,
	(HHT)	changes in frequency. De-	when discriminating com-	89, 90, 91]
		composition does not re-	ponents in narrow band
		quire a window selection.	signals. Time-frequency
					decomposition is more
					difficult to interpret than
					other methods.
	Variational Mode Decom-	Corrects sampling and	Boundary effects are one	[66, 93]
	position (VMD)	noise sensitivity compared	of its main drawbacks.
		to EMD. The modes are
		extracted concurrently.
	Clarke Transform	Power-invariant transfor-	Time domain analysis	[18, 21]
		mation. Converts a three-	only.	No information
		phase system into a two-	provided in the frequency
		phase system, if the sys-	domain.
		tem is balanced.
	Fortescue Transform	Physical and electrical	Time domain analysis	[95,	19,
		interpretability.	Well-	only.	No information	97]
		known tool for analysis of	provided in the frequency
		unbalanced three-phase	domain.
		systems. Low computa-
		tional cost.

Table 2 .

 2 3: Example of common extracted features using different signal processing techniques

	Signal processing	Extracted features
	method
	STFT	Energy, standard deviation, skewness, entropy. Kurtosis and
		maximum value of 50-Hz contour. Entropy and mean value of
		150-Hz contour [99].
	WT	Using db4 as mother wavelet. Minimum and maximum of the
		four-level decomposition. Energy, minimum, maximum, and
		central cumulants of statistical parameters (second, third and
		fourth order) of 50Hz. [100].
	ST	Minimum and maximum values, standard deviation of S-
		contours. Energy and estimated frequency at the maximum
		amplitude from ST matrix [101].

HHT

Energy, entropy, skewness, minimum and maximum of amplitude curve from the first IMF. Standard deviation, skewness and energy of phase curve

[START_REF] Mishra | Advanced signal processing and machine learning techniques for voltage sag causes detection in an electric power system[END_REF]

.

  104], Huang et al. apply Sequence Forward Search (SFS) for various types of PQ disturbance classification, obtaining a reduction in the number of features and a good classification result with low noise sensitivity. Wrappers are on one hand more computationally expensive than filters and less effective for a large number of features. On the other hand, they generally give better results. They should be implemented carefully as they are prone to overfitting.

• Embedded methods incorporate feature selection as part of the model learning process. They combine the advantages of filter and wrapper methods. Some of the most common examples of this category are Decision Trees (DT), Random Forests (RF), Least Absolute Shrinkeage and Selection Operator (LASSO), etc. For instance,

Table 2 .

 2 4: Main properties of classifiers used for power quality classification

	Classifier	Benefits			Drawbacks		References
	Threshold analysis	Rule-based, easy to im-	Generalization capabili-	[18, 90, 107,
		plement.	No learning	ties are not guaranteed.	108]
		needed.			Threshold	values	are
					manually defined.
	Expert System (ES)	Good performances re-	Costly to develop. Re-	[109, 110]
		gardless of the amount of	quires expert knowledge.
		data. Modular structure
		allows to easily add new
		classes.
	Decision tree (DT)	Decision rules obtained	Instability risk.	Small	[105]
		automatically.	Inter-	changes in data can cause
		pretable decision making	large structural changes.
		process.
	Random Forest (RF)	Automatically generate	Difficult to interpret and	[2, 111]
		uncorrelated DT. Glob-	visualize. Computation-
		ally robust and minimal	ally intensive.
		outlier influence.
	Support Vector Machine	Good trade-off between	Data-driven.	Perfor-	[76,	112,
	(SVM)	generalization	perfor-	mance depend on kernel	113]
		mance and complexity.	selection.	Final model
		Fast training.		difficult to interpret.
	Artificial Neural Networks	Good	global	perfor-	Highly data-driven and	[36, 66, 91,
	(ANN, MLP, PNN, ELM)	mances, versatile and	prone to overfitting. Very	114, 115]
		able to represent complex	difficult to interpret.
		non-linear functions.
	Deep Neural Networks	Outstanding accuracy in	Requires	important	[106,	116,
	(DNN, CNN, LSTM)	pattern recognition tasks.	amounts of data. High	117, 118]
		No need of previous fea-	computational cost for
		ture extraction or selec-	training, usually needs
		tion stage. Optimal fea-	specific hardware (GPU).
		ture extraction is learnt
		during training.

Table 2 .

 2 5: Comparative review of Voltage Sag Cause Classification methods in the literature

	Ref.	Data	ac-	Pre-	Feature	Feature	Feature	VS causes	Accuracy Location
		quisition	processing	extraction	selection	analysis
	[90]	Simulation	-	HHT	-	Rule-based	F(1), TE, MS	96%	No
	[91]	Simulation	-	HHT	-	PNN	F(1), TE, MS	98.63%	No
	[76]	Simulation	-	WT	PCA	SVM	F(1), TE, MS,	99.67%	No
								C
	[113]	Field data	RMS	k-means	-	LS-SVM	F(1), TE, MS	92%	No
					SVD
	[2]	Simulation	-	DWT	-	ML-RF	F(4), MS	95%	Yes
	[77]	Simulation	Segm	VQ,DWT	PSO	NBC	F(1), TE, MS	95%	No
	[78]	Simulation	Norm,	L-WT	ALO	SVM, PNN,	F(5), MS	99.6%	No
				Segm			MLP
	[82]	Simulation	Norm	ST	-	ELM	F(1), TE, MS,	100%	No
								C
	[66]	Simulation	TI-W	ST, VMD	-	ELM	F(1), TE, MS	100%	No
	[36]	Field data	ICA	HOS	TLBO	MLP	F(5), TE, MS	99.14%	No
	[106]	Field data	Norm,	-	-	CNN	F(6), TE, MS	97%	No
				2D-image
	[117]	Simulation	Norm	-	-	Bi-LSTM	F(1), TE, MS,	98%	No
		+	Field					C
		data
	[118]	Simulation	Norm	-	-	Bi-LSTM	F(1), TE, MS,	99%	No
		+	Field					C
		data
	[116]	Simulation	Norm,	-	SDADE	CNN	F(1), TE, MS,	99.22%	No
		+	Field	2D-image		(autoen-		C
		data				coder)
	[111]	Simulation	RMS	FLAG,	-	RF	F(4), TE, MS	99.4%	No
					Shaplet
					transform

Table 3 . 1 :

 31 Definition of classes

	Class Voltage sag cause

Table 3 .

 3 2: Real dataset description

	Class Site 1 Site 2 Site 3	Total
	A1	-	35	-	35
	A2	23	103	29	157
	B1	-	-	-	-
	B2	-	-	90	90
	C1	7	-	2	9
	C2	-	-	-	-
	D	93	-	1	94
	Total	123	138	120	385

Table 3 .

 3 3: LV industrial loads

	Load	Rated power	Other characteristics
	Induction motors	22 to 110 kW	Variable inertia and torque
	Variable speed drives	22 to 110 kW	Scalar control
	Three-phase rectifiers	30 kW (DC loads) 6-pulse
	Three-phase loads	30 kW	cosφ = 0.9
	Single-phase loads	5.5 kW	cosφ = 0.6
	Three-phase isolation transformers 100 to 250 kVA	400 V/ 400 V, of type Dyn and YNyn

Table 3 .

 3

			4: Step pattern comparison
		symmetric2	symmetric1	asymmetric
		D (e-02) Z (%) D (e-02) Z (%) D (e-02) Z (%)
	Intra-class	1.41	71.83	1.85	19.79	2.26	40.74
	Inter-class	6.26	90.11	9.48	20.59	6.89	47.15

Table 3 .

 3

	5: Comparison between centroid estimation methods
		D (e-02) Mean distance to centroid
	Euclidean averaging	5.20
	DBA	1.342
		γ = 0.01	0.842
		γ = 0.10	0.859
	soft-DTW	γ = 0.50	0.986
		γ = 1.00	1.108
		γ = 2.00	1.328
		γ = 10.0	4.607

•

  Test sets, 16 in total. Each test set contains 700 synthetic voltage sags with added Gaussian noise with an SNR = [40, 38, 36, ... ,12, 10] dB. For instance, test set no. 1 contains noisy signatures at SNR = 40 dB, test set no. 2 contains signatures at SNR = 38 dB, and so on.

Table 4 .

 4 2: Prediction and confidence index results using synthetic data -Mean distance with Bootstrapping method

	Class	Database size Test set size Synthetic Synthetic	Bootstrapping F1-score (%) NB-KDE (%) RD (%)
	A1	20	80	100	99.99	96.71
	A2	20	80	100	99.82	97.84
	B1	20	80	100	99.99	95.60
	B2	20	80	100	99.97	97.80
	C1	20	80	100	99.99	99.07
	C2	20	80	100	97.54	94.61
	D	20	80	100	98.14	92.26
	Total/Average	140	560	100	99.35	96.27

Table 4 .

 4 3: Prediction and confidence index results using synthetic data -Centroid estimation method

	Class	Database size Test set size Synthetic Synthetic	Centroid estimation F1-score (%) NB-KDE (%) RD (%)
	A1	20	80	100	99.98	97.04
	A2	20	80	100	97.54	96.55
	B1	20	80	100	99.97	94.54
	B2	20	80	100	99.99	98.06
	C1	20	80	100	99.99	99.31
	C2	20	80	100	91.07	94.68
	D	20	80	100	97.29	95.02
	Total/Average	140	560	100	97.97	96.46

Table 4 .

 4 4: Computational cost comparison

	N k	Database size	Bootstrapping Train time (s) Test time (s) Train time (s) Test time (s) Centroid estimation	Test time ratio*
	20	140	-	0.70	1.88	0.03	23.33
	30	210	-	0.98	3.06	0.03	32.66
	50	350	-	1.55	4.21	0.03	51.91
	80	560	-	2.56	6.89	0.03	85.36
	100	700	-	3.08	8.37	0.03	102.81
			*Test time ratio = bootstrapping test time / centroid estimation test time

Table 4 . 5 :

 45 Description of the experiments For experiments No. 1 and 2, the training and testing sets are from the same source (either synthetic or real). In experiments No. 3 and 4 the training and testing sets are from different sources: the reference training databases have 20 signatures per class but the algorithm is tested on the entire dataset of the second data source.The results of experiment 1 are presented in Table4.6. The classification is optimal across the 7 classes, with an F1-score equal to 100% and an average RD index equal to 96.46%.

	No.	Reference database Type Size	Testing set Type Size
	1	Synthetic	140	Synthetic 560
	2	Real	85	Real	300
	3	Synthetic	140	Real	385
	4	Real	85	Synthetic 700

Table 4 .

 4 6: Experiment 1: Training and testing with synthetic data

	Database size Class A1 A2 B1 B2 C1 C2 D	F1-score (%) RD index (%)
	20	A1	80	0	0	0	0	0	0	100	97.04
	20	A2	0	80	0	0	0	0	0	100	96.55
	20	B1	0	0	80	0	0	0	0	100	94.54
	20	B2	0	0	0	80	0	0	0	100	98.06
	20	C1	0	0	0	0	80	0	0	100	99.31
	20	C2	0	0	0	0	0	80	0	100	94.68
	20	D	0	0	0	0	0	0	80	100	95.02
	Average									100	96.46

Table 4 .

 4 7: Experiment 2: Training and testing with real field data

	Database size Class A1 A2 B1 B2 C1 C2 D	F1-score (%) RD index (%)
	20	A1	15	0	0	0	0	0	0	100	89.24
	20	A2	0	133	0	0	4	0	0	84.27	84.26
	-	B1	-	-	-	-	-	-	-	-	-
	20	B2	0	0	0	70	0	0	0	100	98.08
	5	C1	0	0	0	0	4	0	0	66.67	91.78
	-	C2	-	-	-	-	-	-	-	-	-
	20	D	0	0	0	0	0	0	74	100	99.81
	Average									90.19	92.64

Table 4 .

 4 8: Experiment 3: Training with synthetic data, and testing with real field data

	Database size Class A1 A2 B1 B2 C1 C2 D	F1-score (%) RD index (%)
	20	A1	35	0	0	0	0	0	0	97.22	88.43
	20	A2	2	155	0	0	0	0	0	99.35	92.94
	20	B1	-	-	-	-	-	-	-	-	-
	20	B2	0	0	0	90	0	0	0	100	92.80
	20	C1	0	0	0	0	9	0	0	100	86.32
	20	C2	-	-	-	-	-	-	-	-	-
	20	D	0	0	0	0	0	0	94	100	74.63
	Average									99.32	87.02

Table 4 .

 4 9: Experiment 4: Training with real field data, and testing with synthetic data

	Database size Class A1	A2 B1	B2	C1 C2	D	F1-score (%) RD index (%)
	20	A1	100	0	0	0	0	0	0	100	89.80
	20	A2	0	100	0	0	0	0	0	100	93.09
	-	B1	-	-	-	-	-	-	-	-	-
	20	B2	0	0	0	100	0	0	0	100	99.21
	5	C1	0	0	0	0	100	0	0	100	89.90
	-	C2	-	-	-	-	-	-	-	-	-
	20	D	0	0	0	0	0	0	100	100	94.05
	Average									100	93.23

Table 4 .

 4 10: Comparison between two distance-to-class calculation methods and two confidence scores

	Confidence score	Robustness to noise	frequency variations Robustness to fundamental	method Selected
	NB-KDE index	-	--
	RD index	+	++	X
	estimation method Distance-to-class	Low computational cost	Accuracy	method Selected
	with bootstrapping Mean distance estimation	--	++
	with soft-DTW Centroid estimation	++	++	X

3 .

 3 Western Electricity Coordinating Council (WECC-CLM): is a composite model consisting of several models: a static load model, a model for power electronics, and four motor models with different mechanical loads. The composite model also includes a representation of a distribution system (a substation transformer, a shunt reactance, and a line feeder). It was developed by WECC for commercial and residential areas and implemented in Siemens PTI PSS/E, an industry-level simulation software. While the WECC-CLM model provides a very detailed representation, it can be difficult to implement due to the high number of parameters to estimate(131 in total)

Table 5 .

 5 1: Load composition estimation -Load description

	Load	ASD + motor Motor Resistive load Electronic load Lighting
	Reference Power (kW)	22	44	15	15	6.5
	Motor load (pu)	0.6	0.6	-	-	-
	Active Power (kW)	13.4	34.8	15	15	6.5
	Reactive Power (kVAR)	2.6	35.6	0	0.33	-0.25

  • Les capacités de généralisation de ces algorithmes ne sont pas fournies car elles ne sont souvent pas évaluées. Ces méthodes sont entraînées partiellement ou entièrement avec des données obtenues à partir de la même source que les données de test. Cette approche n'est pas toujours réalisable lors du déploiement industriel. D'une part car les données réelles sont difficiles à obtenir en amont, et d'une autre part car les utilisateurs finaux ne sont pas toujours en capacité de ré-entraîner les algorithmes eux-mêmes avec leurs propres données lors du déploiement. Ainsi, la méthodologie doit être généralisable et évolutive pour être appliquée dans différents sites industriels. Idéalement, elle devrait être entièrement entraînée sur des données synthétiques préalablement à son déploiement. • La plupart des indicateurs extraits ou features ne sont pas facilement interprétables du point de vue électrique. Les algorithmes interprétables sont préférables pour des questions d'acceptabilité et de fiabilité, ainsi que pour l'analyse des erreurs de classification. L'algorithme à développer doit ainsi intégrer des connaissances théoriques dans le domaine du génie électrique afin de séléctionner des indicateurs physiques pertinents. • L'extraction d'indicateurs scalaires à partir de signaux non-stationnaires comporte un risque de perte d'information, puisque la dépendance temporelle des formes d'ondes électriques n'est pas prise en compte. Des méthodes alternatives d'extraction d'indicateurs doivent être étudiées afin de prendre en compte cet aspect. • Peu de méthodes proposent une approche combinant classification et localisation. La localisation relative des perturbations est essentielle pour les clients industriels, et il convient de l'aborder également.

1 .

 1 Acquisition de données. Deux sources de données sont à disposition : des données de simulation numérique (synthétiques) et des mesures réelles sur le terrain. Pour la génération des données synthétiques, un modèle de réseau industriel est développé sur le logiciel EMTP. 2. Prétraitement. Cette étape assure l'uniformité des données en termes de longueur et de fréquence d'échantillonnage. Il permet également de filtrer les creux de tension incomplets. 3. Extraction d'indicateurs. Les formes d'onde de tension et de courant sont transformés, et des signatures sous la forme de séries temporelles quadridimensionnelles sont extraites par l'application de la transformée de Fourier à court terme (STFT) et de la transformée de Fortescue. L'interprétabilité électrique de ces signatures est essentielle car elle permet de comprendre la prise de décision de l'algorithme, ce qui constitue un avantage par rapport aux autres méthodes dans la littérature. La description de signatures appartenant aux sept classes définies sont également présentées en détail. 4. Analyse d'indicateurs. Nous avons défini une mesure de distance pour comparer les séries temporelles multivariées. Un alignement spatio-temporel est réalisé au préalable en se servant de l'algorithme Dynamic Time Warping (DTW), comme illustré sur la Figure 3. Puis, une stratégie basée sur la distance entre signatures multivariées est utilisée afin d'effectuer la classification. Il s'agit d'une version améliorée de l'algorithme 1NN-DTW, dont les principaux inconvénients sont la sensibilité aux valeurs aberrantes et un coût calculatoire trop élevé. Nous proposons donc un classificateur de type Nearest Neighborhood qui, au lieu de sélectionner le voisin le plus proche pour la classification d'une nouvelle signature, sélectionne la classe la plus proche. Le calcul de la distance d'une nouvelle signature à une classe particulière peut être obtenu : en estimant la valeur moyenne de la distance à toutes les signatures de la classe (mean distance with bootstrapping), ou en estimant la distance aux centroïdes des classes (centroid estimation with soft-DTW ). Enfin, nous proposons de calculer deux indices de confiance associés à la prédiction.

damentale. La méthode est robuste aux niveaux de bruit jusqu'à un SNR = 15 dB et aux variations de fréquence fondamentale jusqu'à une valeur de ϵ = ±0.5Hz. De plus, un indice de confiance sur la prédiction est proposé, augmentant la fiabilité de l'algorithme. Le système offre une mise en oeuvre facile en milieu industriel sans avoir besoin de données enregistrées au préalable. La méthode présente l'avantage d'utiliser une base de données de référence de taille réduite, entièrement composée de données synthétiques. Les principaux avantages de la méthode proposée sont ses capacités de généralisation et la possibilité de déclencher une alerte basée sur l'indice de confiance. La précision de classification obtenue sur des données synthétiques comportant sept causes est de 100%.

Other events such as switching on of large loads are not included in the list because they cause long-term voltage drops (undervoltages), which are not considered as voltage sags.

The thyristor-based bridge is another possible architecture, however its frequency spectrum is wider, and more variable as it depends on the firing angle.

Frequency spectrum for a rectifier with no power factor correction (PFC).

We distinguish here the term "classification" as the identification or labellisation of an event, from the machine learning notion of "classification" associated with supervised learning.

See glossary

Fuzzy logic is an approach that imitates human-like reasoning in which the true value of a number may be between 0 and 1.

Sampling frequency imposed and set as default on the monitoring device.

A feeder line of "feeder" is a power line that carries electricity from a substation to consumers.

The overlap between two successive windows corresponds to W L -W H .

We have implemented the python library "dtw-python" described in[START_REF] Giorgino | Computing and Visualizing Dynamic Time Warping Alignments in R : The dtw Package[END_REF], with minor adjustments.

The signatures correspond to the synthetic dataset of 700 sags for the seven defined classes.

See Table3.2 for a description of classes.

The predicted label is obtained based on the minimum distance to each class as defined in equation(3.14). Therefore, the F1-score is calculated independently of the confidence index.

In a distance-based approach such as Nearest Neighbors, there is not such a training processes that involves the optimization of internal parameters of the classifier. The reference database is sometimes referenced in the text as training data, but this term is only used to facilitate the reader's understanding.

A confusion matrix is commonly used to evaluate the performance of a classification model. It compares the predicted values of the classifier with the actual ground truth values.

Load composition is used to to estimate the proportion of load types in the sub-network. It is different form load class estimation, which is devoted to determine the type of sub-network: residential, industrial, or commercial.

A partir de cette étude, nous proposons une approche dans laquelle nous estimons d'abord la composition de la charge d'un site industriel lors d'un creux de tension sans impact, avant d'étudier le cas avec un creux de tension causant la déconnexion des charges. Les résultats, même prometteurs, ont montré certaines limites de la méthode, notamment en termes d'interactions harmoniques entre les charges. Finalement, nous préconisons plusieurs propositions pour améliorer la méthode.En conclusion, l'essentiel de ce travail de thèse est consacré à la classification des causes de creux de tension et à leur localisation relative par rapport au point de mesure. Bien qu'il existe une grande variété de solutions permettant d'atteindre une précision élevée (entre 92% et 100%) pour la classification des causes de creux de tension, ces méthodes rencontrent certaines limitations qui empêchent leur application dans un contexte industriel réel, telles que la nécessité d'un accès à de grandes quantités de données d'entraînement, la faible interprétabilité du processus de prise de décision, et des capacités de généralisation faibles ou non évaluées sur différentes sources de données. Notre objectif était de proposer une solution permettant d'atteindre des niveaux de précision équivalents ou supérieurs, tout en surmontant ces limitations.

• Concernant l'estimation de la composition des charges auto-déconnectées suite à un creux de tension 1. Pour surmonter certaines des limites du modèle de charge, d'autres techniques d'apprentissage automatique telles que les ANN ou les DNN pourraient être explorées. En effet, ces méthodes offrent plus de flexibilité pour modéliser les interactions harmoniques entre les charges. Il faudrait également utiliser des modèles qui tiennent compte de la puissance de court-circuit et de l'impédance équivalente du réseau amont. 2. Des techniques d'apprentissage automatique ou d'optimisation stochastique pourraient être mises en oeuvre pour l'étape d'estimation de paramètres. Ces méthodes permettraient d'intégrer les incertitudes des paramètres électriques/mécaniques de charges, dues à la grande diversité d'équipements au sein des catégories de charges considérées. 3. Il conviendrait aussi d'explorer les techniques d'agrégation des charges. Pour ceci, il existe quelques propositions dans la littérature, notamment pour l'agrégation des moteurs à induction [153, 154, 155, 156].4. Finalement, des représentations plus réalistes des courbes de tolérance VTC devraient être étudiées, par exemple, avec des fonctions de densité de probabilité comme proposé dans[START_REF] Milanovic | Probabilistic Assessment of Financial Losses due to Interruptions and Voltage Sags-Part I: The Methodology[END_REF].

Remerciements

Classification algorithm scheme

The algorithm consists of four stages: data acquisition, pre-processing, feature extraction, and feature analysis. We distinguish two modes of operation: a reference database constitution mode, prior to the implementation of the algorithm, and a classification mode, during the algorithm's operation. In the first mode, we build a reference database using the processed labeled events. In the classification mode, the unknown voltage sags, are processed by all of the algorithm's steps including the feature analysis step. The two modes are illustrated in Fig. 3.1. powers at the rated voltage V 0 . V is the variable voltage of the device and the coefficients a p , a q , b p , b q , c p and c q are to be estimated. These coefficients depend on the load's characteristics. In [START_REF] Fan | Identifying ZIP Coefficients of Aggregated Residential Load Model Using AMI Data[END_REF] a ZIP model was used to represent a feeder serving a residential area (aggregated loads).

2. Exponential model: it relates the power of a load with its voltage supply by an exponential equation as presented in equation (5.2). P 0 and Q 0 are the active and reactive powers at the rated voltage V 0 . The coefficients n p and n q are estimated from the electrical behavior of the load. This model has fewer parameters to estimate compared to the ZIP model, but might have a more limited domain of validity. This type of model was implemented in [START_REF] Leinakse | Estimation of Exponential and ZIP Load Model of Aggregated Load with Distributed Generation[END_REF] for the modeling and load aggregation of transformers in a MV distribution system with high penetration of distributed generation. The exponential model displayed a lower and more consistent estimation error (MSE = 0.010) than the ZIP model (MSE = 0.274) for six synthetic events.

3. Frequency dependent model: it represents the relationship between power with both voltage and frequency of the load's supply. It is obtained by applying a frequency-dependent coefficient to a ZIP or exponential model. The frequency-dependent factor is described in equation (5.3), where f 0 is the nominal frequency, f the frequency at the bus voltage and a f the parameter to estimate. A ZIP model incorporating a frequency-dependent term was used in [START_REF] Ge | An Event-Oriented Method for Online Load Modeling Based on Synchrophasor Data[END_REF] for an event-oriented method for online load modeling for the Illinois Insitute of Technology microgrid. The load model was able to capture oscillation and damping information successfully, obtaining better approximation results (NMSE = 0.010) compared to a classic ZIP model (NMSE = 0.215) and an induction motor (IM) model (NMSE = 0.152).

Dynamic models

Dynamic models, particularly required for stability studies, are designed to represent the dynamic behavior of loads. They can capture short-term or longterm dynamics depending on the selected model. The most common models in this category are [START_REF] Arif | Load Modeling-A Review[END_REF]:

1. Induction Motor (IM) model: is derived from the dynamic equivalent electrical circuit of the induction motor. Thus, it is considered a physicsbased model. In [START_REF] Aboul-Seoud | Dynamic modeling of induction motor loads for transient voltage stability studies[END_REF], the authors compared an IM model and a static exponential model for aggregated induction motor load modeling. The authors concluded that although the static model could in some cases reproduce the response of motor loads to voltage disturbances, only the IM model was able to capture the dynamic behavior by taking into account the load inertia on the system stability.

2. Exponential Recovery Load Model (ERL): is used to represent slowly recovering loads. The model is represented as a set of non-linear first-order equations, as in (5.4). In [START_REF] Paidi | Development and Validation of a New Oscillatory Component Load Model For Real-Time Estimation of Dynamic Load Model Parameters[END_REF] an oscillatory component load model based on a static, exponential recovery and damped oscillatory components was proposed for representing industrial, residential and commercial loads. The authors validated the model's performance using real field measurements.

Composite models

Composite models are a combination of static and dynamic models. They are expected to provide more accurate results benefiting from the qualities of both models. Some of the most widely used composite models are:

1. ZIP + IM: is one of the most flexible and commonly used composite models, for its simplicity and easy application [START_REF] Arif | Load Modeling-A Review[END_REF]. It combines the characteristics of a static model (ZIP) representing an aggregation of static loads with the dynamic model of an induction motor (IM). It was implemented for stability analysis in [START_REF] Renmu | Composite Load Modeling via Measurement Approach[END_REF] and for transient stability using measured data to predict unseen data in [START_REF] Han | A Real Application of Measurement-Based Load Modeling in Large-Scale Power Grids and its Validation[END_REF].

We select the load model structure as well as the load estimation method and validate the approach to determine the load composition of the industrial site when voltage sag occurs.

Load composition estimation after a voltage sag with impact

Once the load model and estimation method are validated in the first step, we consider voltage sags with self-disconnected loads. For this, the load model will be modified to include the voltage tolerance curves (VTC) of the equipment in the industrial facility.

Figure 5.2: Two-step approach for load composition estimation after a voltage sag with and without impact

The two-step approach is described in Fig. 5.2. However, we assume the following hypotheses:

• All the industrial loads can be classified into one of the five categories described in section 5.2. This first approach does not consider less frequent loads such as arc furnaces.