*Review* **Islanding Detection Methods for Microgrids: A Comprehensive Review**

**Muhammed Y. Worku 1,\*, Mohamed A. Hassan <sup>1</sup> , Luqman S. Maraaba <sup>2</sup> and Mohammad A. Abido 1,3,4**


**Abstract:** Microgrids that are integrated with distributed energy resources (DERs) provide many benefits, including high power quality, energy efficiency and low carbon emissions, to the power grid. Microgrids are operated either in grid-connected or island modes running on different strategies. However, one of the major technical issues in a microgrid is unintentional islanding, where failure to trip the microgrid may lead to serious consequences in terms of protection, security, voltage and frequency stability, and safety. Therefore, fast and efficient islanding detection is necessary for reliable microgrid operations. This paper provides an overview of microgrid islanding detection methods, which are classified as local and remote. Various detection methods in each class are studied, and the advantages and disadvantages of each method are discussed based on performance evaluation indices such as non-detection zone (NDZ), detection time, error detection ratio, power quality and effectiveness in multiple inverter cases. Recent modifications on islanding methods using signal processing techniques and intelligent classifiers are also discussed. Modified passive methods with signal processing and intelligent classifiers are addressing the drawbacks of passive methods and are getting more attention in the recently published works. This comprehensive review of islanding methods will provide power utilities and researchers a reference and guideline to select the best islanding detection method based on their effectiveness and economic feasibility.

**Keywords:** microgrid; islanding detection; local islanding; remote islanding; signal processing

#### **1. Introduction**

Distributed generation (DG) integrated with energy storage, and both renewable and non-renewable energy resources providing power to local loads, forms a microgrid [1,2]. Microgrids increase the reliability and resiliency of the grid by regulating the voltage in medium and low distribution networks. They also offer several advantages and benefits, including a reduction in CO2 emission, improving energy efficiency, the integration of renewable sources, energy access to remote and developing communities, and a reduction in power transmission losses [3–7].

A microgrid has two modes of operation, namely, grid-connected and island (standalone) modes [8,9]. In grid-connected mode, the microgrid operates in parallel with the main utility, and the main grid is responsible for smooth operation by controlling the voltage and frequency. In this mode, the DG units forming the microgrid are controlled and operated in the current control mode, called grid following. In the island mode, the microgrid is operated as an independent power island, controlling its own voltage and

**Citation:** Worku, M.Y.; Hassan, M.A.; Maraaba, L.S.; Abido, M.A. Islanding Detection Methods for Microgrids: A Comprehensive Review. *Mathematics* **2021**, *9*, 3174. https://doi.org/ 10.3390/math9243174

Academic Editor: Nicu Bizon

Received: 3 November 2021 Accepted: 3 December 2021 Published: 9 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

frequency. The DG units in this mode are controlled and operated in voltage control mode, commonly called grid forming [10,11].

Microgrid islanding occurs when the main grid power is interrupted but, at the same time, the microgrid keeps on injecting power to the network, which can be intentional or unintentional [12,13]. Intentional islanding is a controllable operation mode required for the maintenance of the main utility, whereas unintentional islanding is an uncontrollable operation caused by regular faults such as line tripping, equipment failure, or other uncertainties in the power system [14–16] and may degrade the power quality, overload the system, damage equipment and cause safety hazards [17–20]. Therefore, detecting the islanding condition and effectively disconnecting the microgrid within a specified time interval from the distribution network is a necessity. Moreover, in the islanding condition, the conventional protection devices might not operate, as the DG units cannot provide the sufficient fault current for its operation [21]. The authors in [22,23] investigated the design and control requirements to safely island a microgrid operating either in grid-connected or island modes.

The IEEE and IEC revise and modify the DG interconnection and islanding codes frequently to accommodate the fast growing renewable integration [24]. The consequences of unintentional islanding can be avoided by safely following the provided standards from the IEEE and IEC. The increase in DG integration makes the need to detect unintentional islanding a hot research topic. Researchers have developed different islanding detection methods (IDMs) to address the challenges associated with unintentional islanding. Many IDMs proposed in the literature claim a high reliability and better accuracy compared to each other.

This paper presents a detailed review of the different IDMs proposed in the literature. The IDMs are studied considering their effectiveness, performances, feasibility and operational capabilities. Their advantages and disadvantages are also critically analyzed. The rest of the paper is organized as follows. Sections 2 and 3 describe the islanding detection standards and the performance evaluation criteria of IDMs, respectively. Section 4 presents the detailed classification of IDMs. Section 5 presents discussion and recommendation, whereas Section 6 concludes the paper.

#### **2. Islanding Detection Standards**

Figure 1 [8] shows a distribution network connected with distributed energy resources (DERs) and energy storages. The islanding phenomena shown by the dotted lines occurs when the power supply from the grid is interrupted. Unintentional islanding degrades the power quality, complicates orderly power restoration and endangers the lives of utility personnel.

From Figure 1:

*PPV* represents PV array generated power;

*PBAT* is the charging and discharging power of the battery storage system;

*PGEN* is the power generated from the diesel generator;

*PLOAD* is the power drawn by the load;

*PGRID* is the power exchanged between the main grid and the microgrid;

*PCC* is the point of common coupling;

*CB* is the circuit breaker.

The IEEE and IEC offer standards on how the DG units are operated and controlled with the main grid. IEEE Std. 1547 [25] defines islanding as a condition in which part of the power system becomes isolated from the rest of the network. Islanding detection is one of the major issues when deciding if a DG unit is being synchronized with a grid. Islanded operation requires fast, precise, and cost-effective IDMs, which does not affect the quality of supply. Thus, detecting the islanding condition accurately and timely are the two most important factors to save a distribution network from collapsing. Operating DERs in island mode are not allowed under existing standards such as IEC 62,116, IEEE 1547, IEEE 929-2000 and AS4777.3-2005 [26]. In fact, the islanding condition should be detected and

the microgrid disconnected from the main grid within 2 s, as described in IEEE 1547 [27]. The standards describe in detail the operation of the DG, such as disconnecting the DG unit within 2 s, monitoring the magnitude and direction of power flow, appropriate control of voltage, frequency and power quality.

**Figure 1.** Grid and island operation modes in a DER based microgrid.

Table 1 shows some common standards for islanding detection, voltage and frequency ranges, along with the required detection time.

**Table 1.** Standards for microgrid islanding.


#### **3. Performance Evaluation Criteria of IDMs**

Power systems with a high penetration of inverter-based resources, such as wind, solar and energy storage, in the distribution network have a reduced inertia, making them prone to an increased risk of frequency instability [31–33]. For a small disturbance at the point of common coupling (PCC), conventional methods fail to detect the islanding condition.

IEEE 1547 will be used to assess the performance of IDMs in this paper. The performances of different IDMs are evaluated on whether they can detect islanding timely, effectively and accurately. Non-detection zone (NDZ), detection time (DT), error detection ratio (EDR) and power quality (PQ) are the most popular performance indices used to evaluate IDMs. These indices are described in detail.

#### *3.1. Non-Detection Zone (NDZ)*

The NDZ represents a region of power imbalance between the power generated by the DG units and that dissipated by local loads where the islanding detection method fails [34]. The non-detection zone is the main performance indicator for the implemented IDM and is the main reason IDMs fail to detect islanding. The term "power mismatch space" is used to describe IDMs that are based on monitoring voltage, frequency or phase deviation, whereas IDMs that inject a disturbance are expressed in the "load parameter space". Figure 2 presents an NDZ based on a passive islanding detection method.

**Figure 2.** Non-detection zone for over/under voltage (UOV) and over/under frequency (UOF) passive islanding detection method [35].

#### 3.1.1. Power Mismatch Space

For a microgrid operating in island mode, the power imbalance between the power generated from the DG units and that dissipated by the local loads affects the voltage and frequency at the PCC. If the imbalance is nearly equal to zero (Δ*P* and Δ*Q* close to zero), the variation of voltage and frequency will not be enough to detect islanding when the microgrid disconnected from the grid [36]. The NDZ in the power mismatch space is defined as the power imbalance Δ*P* and Δ*Q*, which cannot cause voltage or frequency to exceed the normal limit to detect islanding and is given as [37]:

$$\left(\frac{V}{V\_{\text{max}}}\right)^2 - 1 \le \frac{\Delta P}{P} \le \left(\frac{V}{V\_{\text{min}}}\right)^2 - 1\tag{1}$$

$$Q(1 - \left(\frac{f}{f\_{\rm min}}\right)^2) \le \frac{\Delta Q}{P} \le Q(1 - \left(\frac{f}{f\_{\rm min}}\right)^2) \tag{2}$$

where, *V* and *P* are the rated voltage and the rated active power, respectively;

*Vmin* and *Vmax* are the minimum and maximum microgrid voltages, respectively; *Q* is the quality factor;

*fmin* and *fmax* are the minimum and maximum frequencies, respectively.

#### 3.1.2. Load Parameter Space

Equation (3) defines the NDZ in parameter space as:

$$F\_1(\mathcal{cf}, \mathbf{K}, \mathbf{Q}) < \Delta \mathbb{C}\_{norm} < F\_2(\mathcal{cf}, \mathbf{K}, \mathbf{Q}) \tag{3}$$

where *cf* is the chopping fraction, *K* is the accelerating gain, and Δ*Cnorm* is the resonate capacitance in the range of NDZ.

#### *3.2. Detection Time (DT)*

The detection time is defined as the time taken from the beginning of microgrid disconnection till the end of the IDM detecting islanding.

$$
\Delta T = \left. T\_{IDM} - T\_{trip} \right| \tag{4}
$$

where Δ*T* is the run-on time, *TIDM* is the moment to detect islanding, and *Ttrip* is the moment microgrid disconnects from the grid.

#### *3.3. Error Detection Ratio (EDR)*

Due to load switching, or other disturbances that affect measurement parameters to exceed normal limits, IDMs might detect false islanding, called error detection [38]. This is defined as:

$$E = \frac{N\_{error}}{N\_{error} + N\_{correct}}\tag{5}$$

where *E* is the error detection ratio, *Nerror* is the times of error detection, and *Ncorrect* is the times of correct detection.

#### *3.4. Power Quality (PQ)*

Maintaining the power quality of the microgrid is an important index while selecting IDMs. IDMs that inject a disturbance to the system distort the power output and deteriorate the power quality.

#### **4. Classification of Islanding Detection Methods**

Islanding detection techniques are mainly classified into local and remote [39–41]. Local islanding techniques are further classified as passive, active and hybrid techniques, based on non-detection zone (NDZ), detection speed, power quality, error detection rate and efficacy in multiple inverter cases. Passive islanding techniques are widely used by utilities because of their low cost and that they do not degrade the power quality. However, these methods have a large NDZ, and setting the threshold setting is a challenge. To overcome the limitations of the passive technique, different signal processing and intelligent classifiers have been used in the literature. Figure 3 [12,42] presents the detail classification of IDMs, and these techniques are discussed in detail in the following sections.

#### *4.1. Local Islanding Detection Techniques*

Local islanding detection techniques measure the system parameters at the DG site for islanding detection. The measured parameters include voltage, frequency, active power, reactive power phase angle, impedance and harmonic distortion. Local islanding detection techniques are classified as passive, active and hybrid techniques and are described as follows.

**Figure 3.** Classification of islanding detection methods.

#### 4.1.1. Passive Islanding Detection Techniques

This method measures the system parameters and compares them with a predetermined threshold value for islanding detection. The measured system parameters at the DG terminal or PCC include voltage, frequency, phase angle and harmonics. The passive islanding detection techniques working principle is depicted in Figure 4. Passive islanding detection techniques are mostly used by power utilities as they are simple, low cost, do not degrade the power quality and have a fast detection speed within 2 s, as recommended by IEEE 1547. However, these methods have a large NDZ, the error detection rate is high and setting the threshold requires special consideration. Some of the popular passive IDMs are described below.

**Figure 4.** Passive islanding detection technique.

#### Harmonic Detection (HD)

The HD method is based on comparing the total harmonic distortion (THD) measured at the PCC and a predefined THD to detect islanding. When the microgrid is operated in grid-connected mode, the PCC voltage is a normal sine wave, and the harmonics generated by the load and the inverter are negligible. However, during islanding mode of operation, the harmonics produced by the inverter will distort the PCC voltage and, hence, islanding will be detected [42–45].

This method is easy to implement and is also effective for multiple DGs connected to the same PCC with a detection time of 45 ms. However, selecting the threshold is difficult since grid disturbance can cause error detection, and it might fail to detect islanding for loads with a large quality factor *Q* and a large NDZ. *Q* is defined in Equation (6) as [45,46]:

$$Q = R\sqrt{\frac{C}{L}}\tag{6}$$

Over/under Voltage and over/under Frequency (OUV/OUF)

This method works by comparing the PCC voltage and frequency with a predefined threshold voltage and frequency to detect islanding. The microgrid will be disconnected from the main grid if the measured voltage and frequency at the PCC exceed the thresholds. The microgrid disconnection from the grid deviates the frequency and voltage at the PCC due to an active and reactive power mismatch between the power generated in the DG units and dissipated in the loads.

$$
\Delta P = P\_{load} - P\_{DG} \tag{7}
$$

$$
\Delta Q = Q\_{load} - Q\_{DG} \tag{8}
$$

In grid-connected operation, the main grid injects Δ*P* and Δ*Q*, and the balance of active and reactive power will be kept. When islanding occurs, to keep the active and reactive power balance, the voltage and frequency will drift until Δ*P=0* and Δ*Q=0*. By detecting the deviation in voltage and frequency, OUV/OUF can detect islanding [47–49]. The advantage of this method is that it does not affect the power quality and the cost is low. The disadvantages are that it is difficult to predict the detection time and it has a relatively

large NDZ, with a detection time between 4 ms to2s[50]. This method is suitable for microgrids with some power imbalance between DG and loads.

#### Rate of Change of Frequency (ROCOF)

When the microgrid is disconnected from the main grid with a power mismatch, the frequency will change. The ROCOF method works by measuring *df/dt* for a few cycles and comparing it with a setting threshold. Islanding will be detected if the measured *df/dt* exceeds the predefined threshold [51–54]. Compared to OUV/OUF, this method has a fast detection time of 24 ms, is more sensitive and highly reliable. However, it has a high error detection rate for systems with high load switching and fluctuation. Hence, ROCOF is best suited for loads with less fluctuation as it cannot discriminate whether the frequency change comes from load changes or by islanding [55]. An extension of ROCOF that considers the dynamic behavior of the load is proposed in [56,57]. The authors incorporated the exponential static load model to incorporate the dynamic behavior of the load and determine the threshold to detect islanding.

#### Rate of Change of Frequency over Power (ROCOFOP)

This method works by measuring *∂f/∂PL*, where *PL* is the load power, to detect whether or not islanding occurs. It has a lower error detection ratio, smaller NDZ and higher reliability than ROCOF. This method has a detection time of 100 ms and works efficiently for microgrids that have a small power imbalance between the DG units and the load [58,59].

#### Rate of Change of Power Output (ROCOP)

This method measures the changes in the DG power output (*dP/dt*) over a few cycles and compares it with a setting threshold to detect islanding. Generally, a loss of the main grid produces load changes, and *dP/dt* measured after the microgrid is islanded is greater than *dP/dt* measured before the microgrid is islanded. This method has a fast detection, with a detection time of between 24 and 26 ms, and the power imbalance between the DG units and the load does not affect the detection speed.

#### Phase Jump Detection (PJD)

The working principle of PJD is to monitor the phase jump between the inverter's terminal voltage and the current for islanding detection. During grid-connected mode, the inverter's current will be synchronized with the voltage at the PCC using a phaselocked loop (PLL) to detect the zero crossing of the voltage. In islanding operation, since PLL works only at the zero crossing of the voltage, the inverter output current remains unchanged. However, the voltage will have a sudden jump due to the load phase angle. Comparing the measured phase difference with a predefined threshold can detect islanding. The advantages of this method are that it has a fast detection speed with a detection time of between 10 to 20 ms, does not affect the power quality, works for multiple inverters and is easy to implement [60,61]. However, it is difficult to choose thresholds for microgrids with frequent load switching.

#### Voltage Unbalance (VU)

A microgrid disconnected from the main grid changes the topology of the network that, in turn, causes a voltage unbalance at the DG output. This voltage unbalance can be used for islanding detection if it exceeds the setting threshold. The voltage unbalance at the time *t* is defined as:

$$V\mathcal{U}l\_t = \frac{NS\_t}{VS\_t} \tag{9}$$

where *NSt* and *PSt* are the negative and positive sequence voltage amplitudes, respectively.

The authors in [62] proposed a variational mode decomposition method to obtain the principal modes from the measured three phase voltage signal and employed the mode singular entropy to determine the index for islanding detection. This method is not sensitive to system disturbances caused by variation in normal loads and has a detection time of about 53 ms [63]. However, the system harmonics affects the extraction of the negative sequence voltage component, making the threshold calculation difficult. This method has better applications for systems with frequent load fluctuations, such as motor starting and capacitor bank switching.

Table 2 compares the different passive islanding detection techniques described, with respect to their performance evaluation, such as NDZ, DT, EDR and power quality.


**Table 2.** Summary of different passive techniques.

#### 4.1.2. Active Islanding Detection Techniques

The performance of active detection methods is based on the perturbation and observation concept. These methods perturb system parameters such as frequency, voltage, currents and harmonics. In the presence of a stiff grid, the amplitude of the variation at the PCC is negligible since the grid parameters are dominant. However, during the islanding phenomenon, injecting a disturbance at the PCC results in a significant variation in the DG parameters. Figure 5 shows the basic working principle of active islanding detection techniques.

**Figure 5.** Active islanding detection technique.

Compared to passive islanding techniques, active techniques have a reduced NDZ and low error detection rate. However, active techniques deteriorate the power quality, and additional power electronic circuits are required to inject the perturbations. To observe the effect of perturbation, additional detection time is required, which can affect the stability of the system. Moreover, most of the active detection methods are developed for inverterbased DG units and are not applicable for synchronous generators. Some of the popular active IDMs are described below.

#### Active Frequency Drift (AFD)

An AFD works by slightly distorting the inverter current waveform injected into the PCC. In grid-connected mode, the voltage and frequency are controlled by the grid and are stable. When islanding occurs, the voltage zero crossing occurs earlier than expected because of the distortion of the injected current waveform. This results in a phase error between the inverter's output current and the voltage, which makes the frequency of the inverter output current drift to eliminate the phase error. This drift in frequency again causes an earlier zero crossing than expected. The frequency of the inverter output current will continue to drift until the voltage frequency measured at the PCC exceeds the threshold of OUF to detect islanding.

An islanding detection method based on a low-frequency current injection disturbance, injected in the conventional *dq* controller of the distribution generator, is proposed in [64]. The frequency deviation is processed using the estimation of a signal parameter via the rotational invariance technique to extract the dominant mode of the oscillations present in the PCC frequency signal to detect islanding. The method has a detection time of 0.12 sec, eliminates NDZ and does not affect the power quality.

The chopping fraction given in Equation (10) describes the distortion of the inverter's injected current [65,66].

$$cf = \frac{2t\_z}{Tvutil} \tag{10}$$

where *tz* is the dead time and *TVutil* is the voltage period.

The advantages of AFD are that it has a small NDZ and is easy to implement, with a detection time within 2 s. The disadvantages are that it degrades the power quality if the injected current is heavily distorted, and the method might also fail to detect islanding for multiple inverter cases.

#### Frequency Jump (FJ)

Frequency jump is a modified version of AFD, as it also inserts dead zones into the current waveform. However, unlike AFD, where dead zones are inserted into every cycle, in FJ, it is inserted in every three cycles. In grid-connected mode, the waveform of the voltage at the PCC is not distorted, despite the inverter's distorted current. During islanding, there will be a variation in voltage frequency that will be used to detect islanding [67]. Similar to AFD, this method might fail to detect islanding for multiple inverters working in parallel.

#### Active Frequency Drift with Positive Feedback (AFDPF)

This method is an extension of AFD and works by applying a positive feedback to increase the chopping fraction, which in turn accelerates the frequency deviation to detect islanding more effectively.

$$cf\_k = cf\_{k-1} + F(\Delta \omega\_k) \tag{11}$$

where *cfk* and *cfk*−<sup>1</sup> are the kth and k − 1th cycle chopping fractions, respectively, and can be positive or negative;

*F* is usually a linear function;

*ω<sup>k</sup>* is the frequency of the kth cycle;

Δ ω*<sup>k</sup>* = ω*k*−<sup>1</sup> − ω0.

As compared to AFD, this method has a small NDZ however, it still affects the power quality [68].

Sandia Frequency Shift (SFS)

SFS is also an extension of AFD and works by applying a perturbation to the frequency of the inverter's voltage with a positive feedback. The modified chopping fraction used in SFS is given in [69] as:

$$
\mathcal{cf} = \mathcal{cf}\_0 + \mathcal{K}(f\_{pcc} - f\_{grid}) \tag{12}
$$

where *cf* <sup>0</sup> is the no deviation in frequency chopping factor;

*K* is the accelerating gain;

*fpcc* is the frequency of the PCC voltage;

*fgrid* is the grid frequency.

In grid-connected mode, the voltage frequency of the PCC is maintained by the grid, even if the method attempts to change it. However, during islanding the chopping fraction increases with the increase of *f* at the PCC, which also increases the frequency of the inverter. The process continues until islanding is detected. This method has a detection time of 0.5 s and has the smallest NDZ compared to other active methods [70,71].

#### Sandia Voltage Shift (SVS)

The working principle of SVS and SFS is similar, in that it perturbs the voltage amplitude of the PCC with a positive feedback to change the inverter's output current and power. In grid-connected mode, the power change does not affect the voltage amplitude of the PCC, whereas in island mode, the power change affects the voltage amplitude, which can be used to detect islanding [72]. SVS is easy to implement; however, its disadvantage are that it lightly degrades the power quality, and the inverter's operation efficiency might be reduced because of the change in the output power.

#### Sliding Mode Frequency Shift (SMS)

SMS perturbs the voltage phase of the PCC with a positive feedback and monitors the frequency deviation to detect islanding. In SMS, the current–voltage phase angle of the inverter is set as [73]:

$$\theta = \theta\_m \sin(\frac{\pi}{2} \frac{f^{k-1} - f\_n}{f\_m - f\_n}) \tag{13}$$

where *θ<sup>m</sup>* is the maximum phase angle at the frequency *fm*, *fn* is the rated frequency, and *f <sup>k</sup>*−<sup>1</sup> is the previous cycle frequency.

In grid-connected mode, the microgrid injects active power to the main grid, and its power factor is close to unity, with the phase angle between the inverter current and the PCC voltage close to zero. During islanding operation, the phase angle of the load and the frequency will vary, and if the frequency variation exceeds the threshold, islanding can be detected. The advantages of the SMS method are that it is easy to implement, is highly effective for multiple inverter systems and has a smaller NDZ with a detection time of about 0.4 s [74]. However, this method reduces the power quality and has an impact on the transient stability of the system.

#### Variation of Active and Reactive Power

This works by varying the injected inverter power and monitors the voltage amplitude and frequency variation for islanding detection. During islanding, the active power generated in the DG units will be dissipated in the loads, and the voltage variation must satisfy Equation (14) to balance the active power between DG and the loads.

$$P\_{DG} = P\_{load} = \frac{V^2}{R} \tag{14}$$

When the voltage variation exceeds the threshold of OUV, islanding can be detected. Similarly, the frequency will be affected by the reactive power disturbance, and islanding can be detected when the frequency variation exceeds the threshold. This method is easy to implement and has a small NDZ with a detection time between 0.3 s and 0.75 s. However, the method greatly affects the power quality and transient stability since it varies the inverter output power continuously. The method also might not work effectively for multiple inverters working in parallel.

#### Negative-Sequence Current Injection

The basic working principle of this method is to perturb the three-phase voltagesourced converter with a negative-sequence current and monitor the negative-sequence voltage at the PCC for islanding detection. During normal grid-connected operation, the injected negative-sequence current will not affect the PCC voltage and will flow into the grid since the grid has low impedance. However, during islanding operation, the injected negative-sequence current will flow into the load, creating an unbalance in the PCC voltage, and islanding can be detected if the voltage unbalance exceeds the threshold. The advantages of this method are that it has a very short detection time of 60 ms (3.5 cycles), it is insensitive to load change, has no NDZ and has a higher accuracy than positive-sequence voltage detection [75,76].

#### Impedance Measurement (IM)

This method works by changing the amplitude of the output inverter current. During islanding operation, the current perturbation varies the voltage, and this variation will be calculated as *dv/di*, an equivalent impedance seen from the inverter. Islanding can be detected if the impedance calculated exceeds the threshold [77]. This method has a detection time of between 0.77 s and 0.95 s and has a small NDZ for a single DG system. However, the detection efficiency will decrease in multiple-inverter cases, and the setting the impedance threshold is difficult since it requires the exact grid impedance.

#### Detection of Impedance at Specific Frequency

This method works by injecting specific frequency harmonics and is a special case of harmonic detection method. During grid-connected mode, the injected harmonic current will not affect the PCC voltage and will flow into the grid. During islanding operation, the injected harmonic current will flow into the local load and produce a harmonic voltage at the PCC. Islanding can be detected if the produced harmonic voltage is large enough. The disadvantage of this method is that it affects the operation of equipment such as transformers.

Table 3 describes the different active islanding techniques found in the literature.


**Table 3.** Summary of different active techniques.

#### 4.1.3. Hybrid Islanding Detection Techniques

Hybrid islanding detection techniques are developed from the combination of passive and active detection techniques, and are implemented with two steps. The first step utilizes a passive technique, primarily to detect islanding. If islanding is suspected in the first step, an active technique is employed to accurately detect the islanding [78–81]. Figure 6 depicts the flow chart of the hybrid islanding detection technique.

**Figure 6.** Hybrid islanding detection technique.

The performance indices will improve from the combination of these methods; they generally have a small NDZ, and the power quality degradation is low. However, the cost of the system is high, and the method is time consuming, which makes their real implementation infeasible. Authors in [82,83] described the recent literature on hybrid islanding detection techniques. Some of the hybrid detection methods discussed in the literature are described below.

#### Voltage Unbalance and Frequency Set-Point

This works by combining the voltage unbalance and the positive feedback-based methods. Computing the average voltage unbalance caused by changes in the system and load is the first step of this method. To differentiate whether the voltage unbalance is caused by islanding or system variation, the second step employs a positive feedback-based method to lower the frequency set point gradually, from 60 to 59 Hz in one second, if the measured voltage unbalance is greater than 35% of the average voltage unbalance [84]. If the nominal frequency is maintained at the PCC, then islanding was not the cause of the voltage variation. However, islanding will be detected if the frequency falls below 59.2 Hz in the following 1.5 s. This method has a detection time of 0.15–0.21 s and works best for microgrids with a low penetration of non-synchronous generation units.

#### Voltage Change and Power Shift

This works by combining the rate of change of voltage and the variation of active power methods. Firstly, to suspect islanding, the average rate of change of voltage is calculated for five cycles. If this calculated voltage exceeds the predefined threshold value, islanding is suspected, and the second method injects extra active power into the system to confirm whether or not islanding has occurred. For normal grid-connected operation, the grid regulates the PCC voltage to be within a predefined interval and compensates the extra injected active power. However, if the microgrid is in islanding operation, the extra active power will increase the voltage amplitude at the PCC. Therefore, islanding can be detected by monitoring the PCC voltage with a detection time of less than 0.5 s. Using reactive power instead of active power for the second step was proposed in [85–87]. The injected reactive power causes an increase in the PCC frequency and islanding can be detected if the frequency variation is more than a predefined threshold.

#### Voltage Fluctuation Injection

This works by combining the rate of change of frequency and voltage methods to detect islanding [88]. The rate of change of voltage and the rate of change of frequency at the PCC are monitored as a first step to detect if one of them exceeds the predefined threshold to indicate islanding might have occurred. Then, the second step employs a voltage perturbation by applying a periodically switching high-impedance load for confirmation. During normal grid-connected operation, the grid stabilizes the PCC voltage perturbation caused by the switching of the high impedance load. However, during the islanding operation, the effect of the periodic perturbation is observed at the PCC voltage to detect islanding. This method does not depend on quality factor and has a detection time of less than 0.216 s, but it might be less efficient for large scale DG units [89,90].

#### Hybrid Sandia Frequency Shift and *Qg*−*f* Method

This technique modifies the Sandia frequency shift method to reduce the NDZ by adding a *Qg*−*f* droop curve to maintain the optimal gain *Kf* at a stable value [91]. The optimal gain *Kf* is directly proportional to the quality factor of the load; however, when the quality factor is more than five, *Kf* will be too large to create a false detection and can even cause system instability. The authors in [88] proposed a *Qg*−*f* droop curve method to keep *Kf* to a safe value and monitor the change of frequency for islanding detection. During normal operation, the reactive power is controlled by the grid; however, during islanding operation, since the DG units operate at unity power factor and produce no reactive power, this creates a frequency difference between the actual and the rated system frequency. This method monitors this change in frequency for islanding detection and has a detection time of 1.4 s.

#### Rate of Change of Reactive Power and Load-Connecting Strategy

This works by combining the change of reactive power and load connection to detect islanding. If the change in reactive power monitored in the first step is more than the predefined threshold, the second step varies the reactive power by connecting an appropriate load to the microgrid [92,93]. During normal grid-connected operation, the grid regulates the reactive power at the PCC, and the rate of change of reactive power is small. However, during islanding operation, any load change affects the generated reactive power from the DG units. The extra connected load affects the rate of change of reactive power and is used to detect islanding. This method can effectively detect islanding, even in the presence of a small load change, and has a detection time of 40 ms. The disadvantage of this method is that selecting the appropriate extra load is not straightforward.

Table 4 presents the different hybrid islanding detection methods.


#### **Table 4.** Summary of different hybrid techniques.

#### *4.2. Remote Methods*

The remote methods utilize advanced signal processing and communication infrastructure for islanding detection. Remote methods do not have a non-detection zone (NDZ), error detection can be eliminated, and they do not affect the power quality; therefore, they are very sound approaches for islanding detection. Whereas remote methods tend to be

expensive to implement for small microgrids, they are very beneficial for large microgrid applications. Some of the remote methods described in the literature are discussed below.

#### 4.2.1. Power Line Carrier Communication (PLCC)

In the PLCC method, transmitters and receivers are set at the grid and DG side, respectively. Transmitters produce a communication signal along with the power line, and if this signal is interrupted, it indicates that islanding has occurred [94]. The PLCC method has a signal period of four consecutive cycles, and the method can detect islanding if the signal is lost in three consecutive periods. This method has a detection time of 200 ms, has no NDZ, has no impact on power quality and is proven to work on multiple inverter system [95]. However, the transmitter set at the grid side is costly, and this method is not economically feasible for low density DG systems.

#### 4.2.2. Signal Produced by Disconnect (SPD)

Similar to the PLCC method, this method also uses signal transmission between the grid and the DG units to detect islanding. However, this method utilizes microwave, telephone and other forms of signal transmission, rather than the power line. This method has no NDZ; however, it needs more investment to set up the communication line. The SPS method can be extended to add additional control of the DG by the main grid, such as coordinating the power generated between the DG units and the main grid, which is beneficial to black start.

#### 4.2.3. Supervisory Control and Data Acquisition (SCADA)

The SCADA system is based on monitoring main grid parameters such as voltages, currents and frequency, which can also be used for monitoring the status of the circuit breakers and sending them to the microgrid. With proper installation, the NDZ can be eliminated with better efficiency. However, the detection speed of this method is slow and requires a high investment to install a separate instrumentation link. Similar to SPD, this method allows additional control of DG by the main grid.

Table 5 describes the different remote islanding detection techniques.


#### **Table 5.** Summary of different remote techniques.

#### *4.3. Passive and Signal Processing*

The local, remote and hybrid IDMs discussed above have their own advantages and disadvantages. However, accuracy, high detection speed and detecting islanding for multiple inverter system are the most unresolved issues that need more research. Compared to local IDMs, remote-based IDMs are highly reliable and have a negligible NDZ. However, they have a high cost, are complex for implementation and are not preferred [96]. Similarly, passive techniques do not perform well in inverter-based DGs with different system configurations. Passive IDMs are better in terms of them not degrading power quality, their fast detection time and being compatible for all DG types. However, they have a large NDZ, and selecting the threshold is not straightforward in most cases. Modifying passive IDMs using advanced signal processing in time-domain, frequency-domain, and time-frequency-domain can address these two limitations.

These modified passive islanding detection methods improve detection time, reduce NDZ and improve detection performance [97–99]. The signal processing tools help extract and analyze the measured signal in order to perform the required operation. However, these methods require a time consuming data training process. The following section describes the signal processing techniques used for islanding detection.

#### 4.3.1. Fourier Transform (FT)-Based Method

This method extracts the features of a signal at specific frequencies using a frequency domain. Short time Fourier transform (STFT), discrete Fourier transform (DFT) and fast Fourier transform (FFT) have evolved from FT to develop efficient and fast IDMs [100–104]. This method has some limitations, such as low-frequency resolution and reduced spectral estimation [105].

#### 4.3.2. Wavelet Transform (WT)-Based Method

To extract the features of a distorted voltage, current or frequency signal, a wavelet transform is the best signal processing method [106]. This method compares the measured signal wavelet coefficients with a predefined threshold value, and islanding will be detected if these values are larger than the predefined values. The disadvantages of this method are that it can only analyze low frequency band, selecting the threshold is not straightforward, and that the different sampling frequencies and the mother wavelet selection have an impact on the wavelet transform. To analyze the high frequency components, the wavelet packet transform (WPT) is applied using the *d-q* axis of three-phase apparent power [107].

#### 4.3.3. S-Transform (ST)-Based Method

This method is an extension of the WT method and converts a time-domain function into a two-dimensional frequency-domain function. The ST method generates the *S*-matrix and the equivalent time-frequency contours from the measured PCC voltage or current signals. To detect islanding, the spectral energy content of the time-frequency contours is calculated containing the frequency and magnitude deviations [108]. However, the ST method requires more computational memory, and the processing time is large compared to other similar techniques.

#### 4.3.4. Time-Time Transform (TTT)-Based Method

This method works by transforming a one-dimensional time-domain signal into a twodimensional time-domain signal by giving a time–time distribution in a particular window. This method distributes the low and high frequency components in different positions. The TTT method performs well for noisy signals, as the method allows a time-local view of the signal through the scaled window [109].

#### 4.3.5. Auto Correlation Function (ACF)-Based Method

This method measures the power or energy signals and extracts information using finite summation limits. The authors in [110] proposed a modified ACF of the modal current envelope to extract the transient features of sample variance that will be used to detect islanding.

#### 4.3.6. Kalman Filter (KF)-Based Methods

This method uses measured voltage or current signals to extract harmonic features using a time–frequency domain. The authors in [111] proposed a voltage harmonics and selected harmonic distortion (SHD) technique, based on KF, for islanding detection. The residual signal and SHD are used as a condition for islanding detection, where the residual signal is used for the islanding detection and the SHD is used for timely detection.

#### *4.4. Intelligent IDMs*

Passive islanding techniques combined with artificial intelligence provide the most effective and economical method to detect islanding. They have a high accuracy, high reliability, are less complex and have a higher computational efficiency than other methods.

Intelligent IDMs do not require threshold selection. Different intelligent IDMs associated with signal processing, such as artificial neural network (ANN), fuzzy logic (FL), support vector machine (SVM), decision tree (DT) and probabilistic neural network (PNN), are commonly used for islanding detection. The only disadvantage of intelligent IDMs is that they suffer from a large computational burden. Figure 7 [12] shows the basic operation of an intelligent IDM. The method starts offline using a training algorithm to train the system from the PCC measured voltage or current signal. Then, to make the final decision, the online process is activated using an intelligent classifier. Some of the intelligent IDMs are discussed below.

**Figure 7.** Intelligent IDMs.

#### 4.4.1. Artificial Neural Network (ANN)-Based Method

The ANN-based approaches extract important features from the measuring data, which are used for identifying variations in power system parameters [112,113]. ANNbased IDMs perform well for multi-inverters and have a high accuracy and efficiency, but the data processing time is large [114]. The authors in [115] reported an islanding detection technique based on an adaptive neuro-fuzzy inference system (ANFIS).

#### 4.4.2. Decision Tree (DT)-Based Method

This method, with a combination of WPT or discrete wavelet transform (DWT), are used for islanding detection in different intelligent IDMs [116]. WPT or DWT are used to extract information from the measured voltage or current signals, and then the DT classifier further processes these features to detect islanding as tested on the CIGRE distribution system [117].

#### 4.4.3. Probabilistic Neural Network (PNN)-Based Method

This method uses artificial neural hardware in traditional pattern recognition schemes. It can compute non-linear decision boundaries based on a Bayesian classification technique and has four layers as the input, pattern, summation and output [118]. These layers do

not need learning and can perform their functions. PNN-based IDMs are reliable for islanding detection.

#### 4.4.4. Support Vector Machine (SVM)-Based Method

The PCC-measured voltage or current signals are used to extract signature features to indicate islanding occurrence using the SVM classifier and autoregressive modelling [119,120]. IDMs based on SVM have a fast detection speed and high accuracy, but the data training and the algorithm make it too complex for practical application.

#### 4.4.5. Fuzzy Logic (FL) Based Method

FL using DT transformation shows a promising and efficient result in islanding detection [121]. However, the disadvantages of FL is that, because of maximum and minimum combinations, the fuzzy classifiers are highly abstract and sensitive to noisy data [122].

#### **5. Discussion and Recommendation**

The IDMs based on different techniques discussed above have many critical technical issues. These issues have to be fixed to improve their performance and make the IDMs more reliable and efficient. IDMs based on a threshold setting have an NDZ that is hard to eliminate, whereas methods based on a disturbance injection degrade the power quality. On the other hand, signal processing-based techniques have a higher precision and are more robust, versatile, reliable and efficient than the existing IDMs. However, they have a high computational burden. There is always a trade off when selecting the performance indices of the IDMs, and it is important to fully consider the practical operation of the microgrid, with all key performance indices taken into account while selecting the appropriate IDM. Researchers should focus on improving the performance of signal processing and intelligent classifier techniques to come up with the best IDM with a high detection speed, smaller NDZ and lower error detection ratio that does not degrade the power quality.

#### **6. Conclusions**

A comprehensive review of various islanding detection methods (IDMs) is presented in this paper. IDMs are broadly classified into two types: remote and local. Remote-type IDMs use communication signals between the microgrid and the main grid and are fast, reliable and effective with zero non-detection zones. These techniques do not degrade the power quality and can be applied to multi-inverter microgrids; however, they are complex and expensive. On the other hand, local methods are classified as passive, active and hybrid. Passive-based IDMs measure microgrid parameters such as voltage, current, frequency and phase angle and monitor their changes to detect islanding. Passive methods are preferred as they are easy and cheap for practical implementation. However, passive techniques have a large non-detection zone. Active-based IDMs inject a perturbation into the system that affects the power quality, whereas hybrid techniques are a combination of passive and active techniques. Active and hybrid techniques need additional devices to introduce the perturbation, which increases the complexity and implementation cost. Compared to the passive, active, and hybrid techniques, IDMs based on signal processing have been gaining more attention recently for islanding detection to detect the islanding condition accurately and precisely within the shortest period without affecting power quality. Moreover, they have the potential of working for multiple inverters and can also overcome the non-detection zone and threshold setting requirements of conventional techniques. Several methods have been studied and a comparison has been provided based on important performance indices including NDZ, detection time, error detection ration and power quality.

**Author Contributions:** M.Y.W. contributed in identifying and classifying the IDMs, manuscript writing and concluding; M.A.H. contributed in identifying the IDMs, manuscript writing and revising. L.S.M.; contributed in classifying IDMS, manuscript writing and editing; and M.A.A. participated in revising the manuscript, editing and concluding. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors acknowledge the support provided by the Interdisciplinary Research Center for Renewable Energy and Power Systems (IRC-REPS), Research Institute, King Fahd University of Petroleum and Minerals, through Project #INRE2111.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Optimising Energy Management in Hybrid Microgrids**

**Javier Bilbao \*, Eugenio Bravo , Olatz García, Carolina Rebollar and Concepción Varela**

Applied Mathematics Department, Engineering School of Bilbao, University of the Basque Country (UPV/EHU), Pl. Ing. Torres Quevedo, 1, 48013 Bilbao, Spain; eugenio.bravo@ehu.eus (E.B.); olatz.garcia@ehu.eus (O.G.); carolina.rebollar@ehu.eus (C.R.); concepcion.varela@ehu.eus (C.V.)

**\*** Correspondence: javier.bilbao@ehu.eus

**Abstract:** This article deals with the optimization of the operation of hybrid microgrids. Both the problem of controlling the management of load sharing between the different generators and energy storage and possible solutions for the integration of the microgrid into the electricity market will be discussed. Solar and wind energy as well as hybrid storage with hydrogen, as renewable sources, will be considered, which allows management of the energy balance on different time scales. The Machine Learning method of Decision Trees, combined with ensemble methods, will also be introduced to study the optimization of microgrids. The conclusions obtained indicate that the development of suitable controllers can facilitate a competitive participation of renewable energies and the integration of microgrids in the electricity system.

**Keywords:** hybrid microgrids; renewable energies; energy management; electricity system

#### **1. Introduction**

In recent years, the microgrid paradigm has emerged, introduced in the early 21st century by Lasseter [1] as an approach that considers generation and associated loads as a subsystem or microgrid. A microgrid can be considered as a set of loads, generators and storage that can be managed in isolation or connected to the rest of the grid in a coordinated manner to supply electricity reliably [2–7]. In emergency situations (faults, disturbances, etc.), the generators and the corresponding loads can be separated from the distribution grid, maintaining service without damaging the integrity of the system. Although originally associated with electricity grids, the concept has been broadened to any set of equipment, such as loads, storage systems and generators, which operates as a unique manageable system that can provide both electrical and thermal power or fuel to a given area [8]. Today, the operation of Distributed Energy Resources (DER) together with manageable loads (domestic consumption or electric vehicles) and various forms of storage such as batteries, supercapacitors or flywheels, is at the core of the hybrid microgrid concept [9]. A microgrid can operate interconnected to the utility through the main distribution grid, using the so-called Point of Common Coupling (PCC), or in island mode, and can also be interconnected with other microgrid systems, which can lead to more complex structures.

The management of hybrid microgrids presents many challenges [10,11], as they can operate either in island mode or connected to the main grid through the PCC. Proper management of the microgrid is, therefore, necessary for stable and economically efficient operation in both situations. The management system must control and adjust both frequency and voltage in either operating mode, share all the loads between the different Distributed Generators (DG) and storage, control the flow with the main grid and optimize operating costs. In grid-connected mode, voltage and frequency will be set by the main grid, which has synchronous generators and large rolling storage systems.

A necessary step in the difficult process of managing a hybrid microgrid is the mathematical modeling of the aspects (power flow, generation, storage, etc.) of that grid, for

**Citation:** Bilbao, J.; Bravo, E.; García, O.; Rebollar, C.; Varela, C. Optimising Energy Management in Hybrid Microgrids. *Mathematics* **2022**, *10*, 214. https://doi.org/10.3390/ math10020214

Academic Editor: Nicu Bizon

Received: 12 December 2021 Accepted: 5 January 2022 Published: 11 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

a next analysis with method. Several methods are, today, applying for this aim: optimizing energy management. Different methods and algorithms have been developed, and are developing [12,13], such as heuristic methods, optimization methods and Machine Learning methods, which will be predominant in the near future since they can compete in accuracy with the others, and because they can be adapted to different topologies, having enough data.

We present in this article a detailed and deep mathematical modeling development, which is not normally found in the bibliography, and is the basis of two heuristic methods (Hysteresis Band Control and control by means of Fuzzy Logic) and the Decision Trees Machine Learning method, combined with ensemble methods, concluding with a comparison of these methods in a microgrid.

This article is organized as follows: the next section presents the main challenges and functions of hybrid microgrids, along with some benefits that can be derived from their correct operation. Section 3 introduces energy management systems and their need to be modeled mathematically. Section 4 then mathematically models the main aspects of a hybrid microgrid to be taken into account when managing them, such as power flow equations, the generation of electrical energy, its storage, power converters, energy consumption and CO2 emissions and other factors. In Section 5, the metaheuristic techniques are developed: they will be defined and their types will be classified according to the process they follow. The main ones (the Hysteresis Band Control and control by means of Fuzzy Logic) will be developed, although at the end of the chapter, mention will be made of some techniques that are either evolutions of other techniques or that are less well known but are widely used due to their good results. Section 6 introduces the Machine Learning method of Decision Trees, their basis, classification and some applications, to continue with the development of Decision Trees. In that section, bagging and boosting ensemble methods are also introduced. Section 7 details the design and experimental results obtained from the comparison operation of a laboratory microgrid. Section 8 discusses some of the open lines of research, and finally ends with the conclusions, where the two main ones are the mathematical modeling compiled in a single article, including hybrid components (renewable energy and storage), and that the Decision Tree method can be applied to the energy management of a hybrid microgrid, but without a great advantage over more classical methods such as Hysteresis Band Control or the application of Fuzzy Logic.

#### **2. Management of Hybrid Microgrids**

The objective of the management and control of a microgrid is to provide the energy demanded by the loads, using generation and storage systems efficiently and reliably, both under regular conditions and when contingencies occur, whether or not there is a connection to the external network.

Hybrid microgrids introduce a number of operational challenges that must be taken into account in the design of their management and protection systems, due to certain particularities that distinguish them from other systems. According to Olivares et al. [14], the most relevant are:


inertia of the system can lead to considerable frequency deviations in the isolated mode of operation.

• Uncertainty. In hybrid microgrids there is greater uncertainty regarding demand and, above all, generation, as the use of renewable energies means that generation is linked to environmental conditions. Therefore, reliable and economic operation must take into account weather forecasting.

Under these circumstances, the management system must ensure reliable operation of the microgrid. The main functions that can be requested from the management system in the microgrid are [14–16]:


#### **3. Energy Management System (EMS)**

The architecture of a system is defined as the fundamental organization of a system, including its components, the relationships between them and the environment and the principles that govern its design and evolution [17,18]. Among the different control architectures, centralized, decentralized and distributed control architectures have been widely used in industry. On the one hand, the centralized implementation stands out for having greater precision, since the control of the process in question is carried out by a single controller, which receives all the signals provided by the network sensors and, after the control process, issues the values to be taken by the different actuators to achieve correct operation [19]. It is therefore a master–slave configuration in which the controller tries to optimize the operation of a set formed by all the subsystems of the network or process, leaving aside the interest that the subsystems themselves may have in optimizing their own operation at the expense of the common good. On the other hand, the decentralized implementation delegates the control problem to several controllers, reducing the computational expense, but also reducing the accuracy of the controller, as input/output data may overlap. Undoubtedly, the great benefit of this type of architecture lies in the ease of implementation with respect to centralized implementation, due to the reduction of a problem into multiple problems of less difficulty. Halfway between centralized and decentralized implementation is the distributed control architecture. In this architecture, there are problems that are related to each other, allowing the coordination of subsystems. This topology is characteristic of microgrids in which each subsystem has a control objective that is different from the others. The higher the number of components in the microgrid, the higher the data traffic between the different controllers and subsystems, thus requiring more bandwidth in the communication system. However, a distributed implementation can reduce data traffic compared to a centralized one, due to the reduced difficulty of the 'local subproblems' that make up the optimization problem.

A fundamental part of a microgrid is the control system, and more specifically, the control strategy or method that will manage the operation of the microgrid in terms of energy generation and demand, so that the energy storage and the external distribution network can satisfy, at all times, the energy balance in the system as a whole. The Energy Management System (EMS) is the system that performs this task, trying to achieve an efficient use of the different components of the microgrid [20,21]. In order to achieve this goal of efficient use, mathematical modeling of the parts of the system is essential.

#### **4. Mathematical Modeling**

#### *4.1. Power Flow Equations*

Each interconnecting component of an electrical network is called a branch or line and links a node *n* to another node *m* in the network. A line can be modeled by its singlephase equivalent circuit. This equivalent circuit accounts for the electrical properties of the conductors (conductivity and insulation) and the physical properties (diameter and distance between conductors). The most commonly used equivalent circuit of a line is the equivalent Π, although there are other models such as *T*. The complex admittance (inverse of impedance) values are represented by the letter *Y* together with an arbitrary number and each bus can have a generator connected to it at a voltage represented by the letter *V* + bus number.

Expressing the magnitudes in complex form, assuming a permanent and balanced sinusoidal regime, the system can be represented, in compact matrix notation, as follows:

$$\mathbf{I} = \mathbf{Y}\_{BUS}\mathbf{V} \tag{1}$$

where **I** is the column vector of currents at each node, **V** is the column vector of voltages at each node and **Y***BUS* is the admittance matrix. The admittance matrix is composed of complex numbers and has well-known properties: it is symmetrical, each element *Ynn* of the diagonal is the sum of the admittances of the equivalent circuits Π that are connected to node *n*, and the off-diagonal elements *Ynm* are the negative of the admittance of the equivalent Π connected between nodes *n* and *m*. Therefore, the admittance matrix is a square matrix of the same dimension as the number of buses. For each current *n* of the column vector, the power of bus *n* can be calculated as a factor of one, as follows:

$$s\_{\rm li} = \frac{V\_{\rm n} I\_{\rm n}^\*}{S\_{\rm base}} = \frac{V\_{\rm n}}{S\_{\rm base}} \left( \sum\_{m=1}^N Y\_{\rm nm} V\_m \right)^\* = \frac{V\_{\rm n}}{S\_{\rm base}} \sum\_{m=1}^N Y\_{\rm nm}^\* V\_m^\* = \sum\_{m=1}^N v\_{\rm n} v\_{\rm m} e^{j\theta\_{\rm nm}} (G\_{\rm nm} - jB\_{\rm nm}), n = 1, \dots, N \tag{2}$$

where *vn* is the modulus of *Vn* in per unit, *θnm* is the angle difference *θn*–*θm*, and *Ynm* is the element *nm* of the admittance matrix *Gnm* + *jBnm* also in per unit (pu). With Euler's formula, one can write the above equation in rectangular coordinates in the complex plane as shown below:

$$s\_{ll} = \sum\_{m=1}^{N} v\_{ll} v\_{m} (\cos \theta\_{nm} + j \sin \theta\_{nm}) (G\_{nm} - jB\_{nm}) \; n = 1, \ldots, N \tag{3}$$

Remembering that the complex part of the apparent power is the reactive power and the real part is the active power, the two can be separated as follows:

$$\begin{aligned} p\_{\text{\textquotedblleft}} &= \sum\_{m=1}^{N} \upsilon\_{\text{\textquotedblleft}m} (\text{G}\_{\text{nm}} \cos \theta\_{nm} + \text{B}\_{nm} \sin \theta\_{nm}) \; n = 1, \dots, N\\ q\_{\text{\textquotedblleft}m} &= \sum\_{m=1}^{N} \upsilon\_{\text{\textquotedblleft}m} (\text{G}\_{\text{nm}} \sin \theta\_{nm} - \text{B}\_{nm} \cos \theta\_{nm}) \; n = 1, \dots, N \end{aligned} \tag{4}$$

The above representation is compact and allows observation of the asymmetric and non-linear character of the power flow equations, but to apply the relevant approximations and obtain a linearization, the special structure of the admittance matrix is considered: the elements of the diagonal *Ynn* are the negative of the sum of the off-diagonal elements (negative of the admittance of the equivalent Π connected between nodes *n* and *m*) of the

corresponding rows and of the shunt admittances of the bus (superscript *sh*). This can be seen in the three bus example below:

$$Y\_{BUS} = \begin{pmatrix} Y\_1^{\text{sl}} + Y\_2 + Y\_8 + Y\_9^{\text{sl}} & -Y\_2 & -Y\_8 \\ -Y\_2 & Y\_2 + Y\_3^{\text{sl}} + Y\_5 + Y\_4^{\text{sl}} & -Y\_5 \\ -Y\_8 & -Y\_5 & Y\_8 + Y\_7^{\text{sl}} + Y\_5 + Y\_6^{\text{sl}} \end{pmatrix} = \begin{pmatrix} Y\_{11} & Y\_{12} & Y\_{13} \\ Y\_{21} & Y\_{22} & Y\_{23} \\ Y\_{31} & Y\_{32} & Y\_{33} \end{pmatrix} \tag{5}$$

In the notation currently considered, with *ynm* = *ymn*, the admittance matrix is as follows:

$$Y\_{\rm BUS} = \begin{pmatrix} y\_{12}^{\rm sh} + y\_{12} + y\_{13} + y\_{13}^{\rm sh} & -y\_{12} & -y\_{13} \\ -y\_{21} & y\_{21}^{\rm sh} + y\_{21} + y\_{23} + y\_{23}^{\rm sh} & -y\_{23} \\ -y\_{31} & -y\_{32} & y\_{31}^{\rm sh} + y\_{31} + y\_{32} + y\_{32}^{\rm sh} \end{pmatrix} \tag{6}$$

$$Y\_{nm} = \begin{cases} -y\_{nm} & \text{if } m \neq n \\ \sum\_{m=1, m \neq n}^{N} y\_{nm} + y\_{nm}^{\text{sh}} & \text{if } m = n \end{cases} \tag{7}$$

$$\sum\_{m=1, m \neq n}^{N} y\_{nm} + y\_{nm}^{\text{sh}} = y\_n^{\text{sh}} + \sum\_{m=1, m \neq n}^{N} g\_{nm} + jb\_{nm} = j \left( b\_n^{\text{sh}} + \sum\_{m=1, m \neq n}^{N} b\_{nm} \right) + g\_n^{\text{sh}} \sum\_{m=1, m \neq n}^{N} g\_{nm} \tag{8}$$

Then, the active and reactive power equations for each node can be rewritten, based on the admittance of each line between bus n and bus m, *Ynm* = −*ynm* = −*gnm* − *jbnm*, and of the shunt of bus *n* as:

$$p\_n = \left(g\_n^{\text{sh}} + \sum\_{m=1, m \neq n}^{N} g\_{nm}\right) v\_n^2 - \sum\_{m=1, m \neq n}^{N} v\_n v\_m \left(g\_{nm} \cos \theta\_{nm} + b\_{nm} \sin \theta\_{nm}\right) n = 1, \dots, N \tag{9}$$

$$q\_{\rm ll} = -\left(b\_{\rm n}^{sh} + \sum\_{\substack{m=1, m \neq n}}^{N} b\_{\rm nm}\right) v\_{\rm n}^2 - \sum\_{m=1, m \neq n}^{N} v\_{\rm n} v\_{\rm m} \left(g\_{\rm nm} \sin \theta\_{\rm nm} - b\_{\rm nm} \cos \theta\_{\rm nm}\right) \, n = 1, \dots, N \tag{10}$$

You can group the summation terms, since the sum is over the same set, and take out common factor *vn* in both equations by grouping the conductance and susceptance coefficients:

$$p\_n = g\_n^{\text{sh}} v\_n^2 + \sum\_{m=1, m \neq n}^{N} g\_{nm} v\_n (v\_n - v\_m \cos \theta\_{nm}) - b\_{nm} v\_n v\_m \sin \theta\_{nm} \, n = 1, \dots, N \tag{11}$$

$$q\_n = -b\_n^{sh} v\_n^2 + \sum\_{m=1, m \neq n}^{N} -b\_{nm} v\_n (v\_n - v\_m \cos \theta\_{nm}) - g\_{nm} v\_n v\_m \sin \theta\_{nm} \, n = 1, \dots, N \tag{12}$$

Different assumptions can now be made to linearize the power flow. Each approach assumes a different approach to the problem; however, the following assumptions are common according to the normal operation of an electrical power system [22,23]:

1. The voltage values, expressed in per unit (pu), are very close to l.

2. The difference between the angles of two interconnected buses is a small number close to 0.

Based on the above, to eliminate the non-linearity of the trigonometric functions, they are approximated by their Taylor series centered at zero and neglecting terms of order equal to or greater than three:

$$
\cos \theta\_{nm} \sim 1 + \frac{\theta\_{nm}^2}{2}, \sin \theta\_{nm} \sim \theta\_{nm} \tag{13}
$$

However, the quadratic term is non-linear and, in addition, there are products of several variables, which is also a non-linear function. The power equations at each bus *n* take the following form:

$$p\_n = g\_n^{sh} v\_n^2 + \sum\_{m=1, m \neq n}^{N} g\_{nm} v\_n \left( v\_n - v\_m - v\_m \frac{\theta\_{nm}^2}{2} \right) - b\_{nm} v\_n v\_m \theta\_{nm} \ n = 1, \dots, N \tag{14}$$

$$q\_n = -b\_n^{sh} v\_n^2 + \sum\_{m=1, m \neq n}^{N} -b\_{nm} v\_n \left(v\_n - v\_m - v\_m \frac{\theta\_{nm}^2}{2} \right) - g\_{nm} v\_n v\_m \theta\_{nm} n = 1, \dots, N \tag{15}$$

At this point, there is currently no consensus on the best approximation of the nonlinear terms. Several approaches consider a second-order approximation based on the Taylor series of the products of two variables, which results in a linear form of the equations, except for the losses. In fact, power losses are non-convex quadratic functions which forces relaxations, such as piecewise linear function approximation or binary expansion, which give rise to a mixed integer linear programming problem (known by its acronym as MILP or MIP), or also convex relaxation, which generates a second-order conic programming problem; both are linear cases considering *v*<sup>2</sup> *<sup>n</sup>* as a variable [24–27].

Since MIP-type problems are more costly to solve in terms of resources and time than continuous linear problems, other authors propose not considering the second order terms of the Taylor series expansion, both of the trigonometric functions and of the products of the variables, in such a way that the resulting flow is symmetric and allows calculating the voltage and angle value at each bus, as well as the powers, however in this approximation, losses are not represented [28].

Taking the first term of the trigonometric functions expansion, and according to Yang et al. [29], three linear approximations of the term *vn*(*vn* − *vm*) are compared, being the one with the least error, in terms of voltage and active power flow, the one that considers a decomposition of the bus voltages as *vn* =1+ Δ*vn*, where Δ*vn* is an order of magnitude smaller than *vn*, therefore:

$$\upsilon\_{\boldsymbol{\eta}}(\upsilon\_{\boldsymbol{\eta}}-\upsilon\_{\boldsymbol{m}}) = (1+\Delta\upsilon\_{\boldsymbol{\eta}})(\Delta\upsilon\_{\boldsymbol{\eta}}-\Delta\upsilon\_{\boldsymbol{m}}) = \Delta\upsilon\_{\boldsymbol{\eta}} - \Delta\upsilon\_{\boldsymbol{m}} + \Delta\upsilon\_{\boldsymbol{m}}\Delta\upsilon\_{\boldsymbol{\eta}} - \Delta\upsilon\_{\boldsymbol{m}}\Delta\upsilon\_{\boldsymbol{\eta}} \sim \Delta\upsilon\_{\boldsymbol{\eta}} - \Delta\upsilon\_{\boldsymbol{m}} \tag{16}$$

In the resulting product of the expansion, the multiplication of the difference of the voltages, Δ*vn*Δ*vn*, is neglected with respect to its nominal value of 1 pu, since the result is at most two orders of magnitude smaller than *vn*; thus *vn*(*vn* − *vm*) *vn* − *vm*. Furthermore, the voltage squared multiplying the shunt terms is simply approximated by the value of the voltage at that bus, and in the case of reactive power the shunt conductance is assumed to be negligible compared to the shunt susceptance. This gives rise to a linear problem in the voltage and angle variables [30,31].

Expanding the product terms of variables by their first order Taylor series gives the following:

$$\begin{array}{l} \upsilon\_{n}\upsilon\_{m}\theta\_{nm} \simeq \upsilon\_{n,0}\upsilon\_{m,0}\theta\_{nm} + (\upsilon\_{n}\upsilon\_{m} - \upsilon\_{n,0}\upsilon\_{m,0})\theta\_{nm,0} = \theta\_{nm} \\ \upsilon\_{n}\upsilon\_{m}\theta\_{nm}^{2} \simeq \upsilon\_{n,0}\upsilon\_{m,0}\theta\_{nm}^{2} + (\upsilon\_{n}\upsilon\_{m} - \upsilon\_{n,0}\upsilon\_{m,0})\theta\_{nm,0}^{2} = \theta\_{nm}^{2} \end{array} \tag{17}$$

where the subscript 0 denotes the point around which the expansion is performed, which in the case of the voltages is *vn*,0 = *vm*,0 = 1, and in the case of the angle difference is *θnm*,0 = *θ*<sup>2</sup> *nm*,0 = 0, justified by the usual operating conditions in most power systems [23]. The power equations at each bus are then as follows

$$p\_n = g\_n^{sh} v\_n^2 + \sum\_{m=1, m \neq n}^{N} g\_{nm} \left( v\_n^2 - v\_n v\_m - \frac{\theta\_{nm}^2}{2} \right) - b\_{nm} \theta\_{nm} \tag{18}$$

$$q\_{n} = -b\_{n}^{sh}v\_{n}^{2} + \sum\_{m=1, m \neq n}^{N} -b\_{nm} \left(v\_{n}^{2} - v\_{n}v\_{m} - \frac{\theta\_{nm}^{2}}{2}\right) - g\_{nm}\theta\_{nm} \tag{19}$$

Neglecting *vnvm* and *θ*<sup>2</sup> *nm*, together with the approximation *v*<sup>2</sup> *<sup>n</sup> vn*, the model of (29) is obtained, while taking the voltage squared as independent variable, the above equations can be rewritten, considering the following voltage product transformation:

$$v\_n v\_m = \frac{v\_n^2 + v\_m^2}{2} - \frac{\left(v\_n - v\_m\right)^2}{2} \tag{20}$$

As

$$p\_n = g\_n^{sh} v\_n^2 + \sum\_{m=1, m \neq n}^{N} g\_{nm} \left( \frac{v\_n^2 - v\_m^2}{2} + \frac{(v\_n - v\_m)^2}{2} - \frac{\theta\_{nm}^2}{2} \right) - b\_{nm} \theta\_{nm} \tag{21}$$

$$q\_n = -b\_n^{sh} v\_n^2 + \sum\_{m=1, m \neq n}^{N} -b\_{nm} \left( \frac{v\_n^2 - v\_m^2}{2} + \frac{(v\_n - v\_m)^2}{2} - \frac{\theta\_{nm}^2}{2} \right) - g\_{nm} \theta\_{nm} \tag{22}$$

Approximating the terms with *θ*<sup>2</sup> *nm* and (*vn* − *vm*) 2 , which account for the losses, by linear functions of *θnm* and *v*<sup>2</sup> *<sup>n</sup>* − *<sup>v</sup>*<sup>2</sup> *<sup>m</sup>*, we would now yield the method proposed by Yang et al. [32], which is linear with respect to the voltage squared and the angle value. However, rearranging the above expressions,

$$p\_n = g\_n^{sh} v\_n^2 + \sum\_{m=1, m \neq n}^{N} g\_{nm} \frac{v\_n^2 - v\_m^2}{2} - b\_{nm} \theta\_{nm} + g\_{nm} \left(\frac{(v\_n - v\_m)^2}{2} - \frac{\theta\_{nm}^2}{2}\right) \tag{23}$$

$$q\_{n} = -b\_{n}^{\text{sh}}v\_{n}^{2} + \sum\_{m=1, m \neq n}^{N} -b\_{nm}\frac{v\_{n}^{2} - v\_{m}^{2}}{2} - g\_{nm}\theta\_{nm} - b\_{nm}\left(\frac{(v\_{n} - v\_{m})^{2}}{2} - \frac{\theta\_{nm}^{2}}{2}\right) \tag{24}$$

and considering the following approximation of(*vn* − *vm*) 2 /2 around the point *vn* = *vm* = 1

$$\frac{\left(v\_n - v\_m\right)^2}{2} \simeq \frac{1}{2} \left[ (v\_n - v\_m) \frac{v\_n + v\_m}{2} \right]^2 = \frac{\left(v\_n^2 - v\_m^2\right)^2}{8} \tag{25}$$

the formulation of (26) is obtained:

$$p\_n = g\_n^{sh} v\_n^2 + \sum\_{m=1, m \neq n}^{N} g\_{nm} \frac{v\_n^2 - v\_m^2}{2} - b\_{nm} \theta\_{nm} + g\_{nm} \left(\frac{\left(v\_n^2 - v\_m^2\right)^2}{8} - \frac{\theta\_{nm}^2}{2}\right) \tag{26}$$

$$q\_n = -b\_n^{sh} v\_n^2 + \sum\_{m=1, m \neq n}^{N} -b\_{nm} \frac{v\_n^2 - v\_m^2}{2} - g\_{nm} \theta\_{nm} - b\_{nm} \left(\frac{\left(v\_n^2 - v\_m^2\right)^2}{8} - \frac{\theta\_{nm}^2}{2}\right) \tag{27}$$

In this case, given that the equations show the value of the angle difference *θnm* and its square *θ*<sup>2</sup> *nm*, as well as the square of the squared voltage difference (*v*<sup>2</sup> *<sup>n</sup>* − *<sup>v</sup>*<sup>2</sup> *m*) 2 , it is necessary to linearize these two terms by piecewise linear functions to properly estimate the losses.

At this point it is necessary to make a decision on the approximation of the power flow equations to be implemented in solving the problem. However, it requires two additional linearizations by piecewise functions, which makes it slower to solve than the simplified method of (30), which does not include terms that account for losses but is computationally faster since it only involves one linearization corresponding to the power limitation of the lines.

Since the objective of implementing these equations is to ensure the correct operation of the power system under study under different scenarios, the procedure involves the resolution of these flows in all cases, as the number of scenarios is very high, a more elaborate description results in possibly unaffordable computation times, which depends on the number of scenarios considered and the time intervals contained in each one.

In order not to limit the scenarios to be considered and the time steps, we implement the model evaluated by Yang et al. [29] and Morvaj et al. [33], originating from (30). As the case considered is a distribution network whose lines have a length of less than 5 km, without committing significant error, the shunt admittance of the lines can be considered null [34,35]. Naming *vnm* = *vn* − *vm*, the linearized power flow equations are as shown below:

$$p\_n = p\_n^G - p\_n^D = \sum\_{m=1, m \neq n}^{N} g\_{nm} v\_{nm} - b\_{nm} \theta\_{nm} = \sum\_{m=1, m \neq n}^{N} p\_{nm} n = 1, \dots, N \tag{28}$$

$$q\_n = q\_n^G - q\_n^D = \sum\_{m=1, m \neq n}^{N} -b\_{nm}v\_{nm} - g\_{nm}\theta\_{nm} = \sum\_{m=1, m \neq n}^{N} q\_{nm}n = 1, \dots, N \tag{29}$$

where the superscripts denote power generation, *G*, and power demand, *D*. Each summand *nm* of the active and reactive power expressions is the power flow per each line linking bus *n* to bus *m*.

On the other hand, the maximum power constraint through the lines is defined by the inequality:

$$p\_{nm}^2 + q\_{nm}^2 \le s\_{\text{max},nm}^2 = i\_{\text{max},nm}^2 \\ v\_{\text{x}}^2 \forall \text{x} = n \text{ }\_{\text{}}m \tag{30}$$

where *Smax*, is the maximum apparent power in per unit (pu) that can circulate through each line, *imax*, is the maximum current that can pass through the conductor and *vx* is the voltage at bus *n* or *m*. The manufacturers provide the maximum current limit because the limiting factor is the temperature of the conductor due to the heat caused by the passage of the current. The model of the conductors is given by the standards and technical requirements that can be found in the technical literature [36], and the value of the maximum current in [37,38], which is 285 A for underground conductors (20 ◦C ground temperature and 70 ◦C conductor temperature), and 262 A for overhead conductors (at 75 ◦C conductor temperature and 35 ◦C ambient temperature). This is what is known as the thermal limit of the conductors and is a factor to be considered in lines shorter than 80 km, as is the case; while between 80 and 320 km the limiting factor is the voltage drop, in lines longer than 320 km it is the stability of the angle [34]. The current limit depends on the conductor temperature, since the electrical parameters of the conductor vary with temperature, however, the change in these parameters would modify the power flow and therefore an iterative calculation would be necessary, in fact, this maximum current decreases with ambient temperature and the dynamics of heat transfer would have to be considered, which gives rise to a non-linear problem: the current limit of the conductors increases non-linearly with conductor temperature and decreases with ambient temperature [38]. This variation is not considered to be significant and the maximum current is assumed to be fixed and equal to that given above, at the given temperature.

#### *4.2. Generation of Electrical Energy*

The parameter representing the photovoltaic active power leaving the solar field and entering the inverter at each bus, in each scenario and time instant, is given by the following formula [39]:

$$P\_{pv}(t) = P\_{nom} \frac{G(t)}{G\_{\text{fl}}} \left[ 1 - a \left( T(t) + \frac{G(t)}{800} [\text{NOCT} - 20] - 25 \right) \right] \tag{31}$$

where *T* is the ambient temperature, *Pnom* is the power at nominal conditions, *Gn* is the nominal irradiance in W/m<sup>2</sup> (this value is sometimes generalized to 1000 [40]), *G* is the incident irradiance in W/m2, α is the parameter of the power-temperature characteristic in %/◦C, *NOCT* takes the value of 45 ◦C and is the nominal operating temperature of the cell at 800 W/m2 with 20 ◦C ambient temperature and 1 m/s wind. The alpha factor of the efficiency is negative, which implies an increase in efficiency with decreasing temperature. The inverter efficiency and other characteristics are introduced in more detail in later sections.

In the case of wind turbines, the parameter denoting the maximum energy extractable from the device is calculated using its power curve, but applying the inverter efficiency to the proportional part of the power flowing through the inverter:

$$P\_{\rm wt}^{\rm max}(t) = P\_{\rm wt}(t)(0.7 + 0.3\eta\_{b2b}(t)) \simeq P\_{\rm wt}(t)(0.7 + 0.3 \cdot 0.965) = 0.9895 P\_{\rm wt}(t) \tag{32}$$

$$P\_{wt}(t) = \begin{cases} 2100 \left[ 1 - e^{-\left(\frac{c(t)}{5.692}\right)^{3.98}} \right] & 2 \text{ m/s} \le c(t) < 20.5 \text{ m/s} \\\\ 4378 - 111c(t) & 20.5 \text{ m/s} \le c(t) \le 25 \text{ m/s} \\\\ 0 & \text{any other case} \end{cases} \tag{33}$$

This gives the upper limits of active power injection to the grid, those of reactive power are related by means of the equations of the converters presented in the section corresponding to power converters.

#### *4.3. Electrical Energy Storage*

This section presents the equations that give the state of charge of each storage system on each bus, in each scenario and at each moment in time. The active and reactive power variables that appear represent the energy output/input of the system through the converter that connects the *ESS* (Energy Storage System) to the grid.

The following stochastic variables are defined for the load power, the discharge power and the state of charge of each *ESS* at each bus *n* for time *t*:

$$P\_{\rm in}^{ESS,n}(t) , P\_{\rm out}^{ESS,n}(t) , \text{SOC}(t) \ge 0 \tag{34}$$

These powers are the real effective powers received or delivered by the *ESS*s, that is, they include the load and converter efficiency in their definition, which is explained in the following development.

The weight of the problem lies in the binary part, not in the linear part, therefore it is convenient to create a variable for loading and another for unloading and to use only a binary variable and its complementary to avoid simultaneous loading and unloading. This restriction is as follows:

$$P\_{\rm in}^{ESS,n}(t) \le b^{ESS,n}(t) P\_{\rm in,max}^{ESS,n} P\_{\rm out}^{ESS,n}(t) \le (1 - b^{ESS,n}(t)) P\_{\rm out,max}^{ESS,n} \,\forall ESS, n \tag{35}$$

where the subscript max denotes the nominal power of each storage system (*ESS*) on each bus *n* and the binary variable is *bESS,n* which takes the value 1 or 0 at each time t for each *ESS* on each bus *n* and each scenario. This constraint forces the load to zero if there is unloading and vice versa.

On the other hand, the load state model is estimated to be linear and without capacity reduction due to unloading depth or gradation. The correctness of this assumption is ensured by adding a number of constraints to limit the state of charge to a safe range of each *ESS*. The state of charge of each *ESS* on each bus n at each time t and scenario is given by the following expression:

$$SOC\_{ESS,n}(t) = (1 - a\Delta t)SOC\_{ESS,n}(t-1) + P\_{\text{in}}^{ESS,n}(t)\Delta t - P\_{\text{out}}^{ESS,n}(t)\Delta t \,\forall ESS, n \tag{36}$$

where *SOC* represents the state of charge, the parameter α is the relative self-discharge per unit time and Δ*t* is the time step. The upper limit of the state of charge of each *ESS* is defined by the following constraints:

$$SOC\_{ESS,\mu}^{\text{max}} DOD\_{ESS,\mu} \le SOC\_{ESS,\mu}(t) \le SOC\_{ESS,\mu}^{\text{max}} (1 - DOD\_{ESS,\mu}) \tag{37}$$

where *DOD* is the relative depth of discharge parameter of each *ESS* on each bus, set to zero for flow batteries and 0.1 for lithium-ion batteries and hydrogen cells, as already introduced in the section on system components.

All these quantities are expressed in International System units, but to link these subsystems to the grid they must be expressed in per unit, and for this purpose they are simply redefined through the quotient between the base power. In addition, the definition of the charging and discharging power variables relate them to the converter that connects these systems to the grid and is where the charging and discharging efficiencies are applied, as well as those of each converter, as noted at the beginning:

$$p\_{\text{out}}^{\text{ESS,n}} = \frac{P\_{\text{out}}^{\text{ESS,n}}}{S\_{\text{base}}} = \frac{p\_{\text{ESS2net,n}}(t)}{\eta\_{\text{out}}^{\text{ESS,n}} \eta\_{\text{conv}}^{\text{ESS,n}}}, \\ p\_{\text{in}}^{\text{ESS,n}}(t) = \frac{P\_{\text{in}}^{\text{ESS,n}}(t)}{S\_{\text{base}}} = p\_{\text{net2ESS,n}}(t) \eta\_{\text{in}}^{\text{ESS,n}} \eta\_{\text{conv}}^{\text{ESS,n}} \tag{38}$$

where *conv* refers to the power converters and the subscripts *ESS*2*net* indicate power transfer from the storage system to the grid and *net*2*ESS* from the grid to the storage system.

In the application case, the charging and discharging efficiencies are those justified in the system components section and the efficiency of the converters is approximated in a constant way based on the commercial model also chosen in the aforementioned section, with a value of 0.98.

Finally, constraints are necessary to force the initial SOC equal to the final SOC, during the time horizon, in order to increase the lifetime of the *ESS*. This variable, *SOC*<sup>0</sup> *ESS*,*n*, is not stochastic:

$$\text{SOC}\_{ESS,n}(t\_0) = \text{SOC}\_{ESS,n'}^0 \text{ SOC}\_{ESS,n}(T) = \text{SOC}\_{ESS,n}^0 \tag{39}$$

#### *4.4. Power Converters*

One aspect to consider for power converters is their performance curve: PV plants and storage systems need to be able to deliver more than a certain percentage of the inverter's rated power for the power output to be effective and work above the bend of the performance curve, whenever it is considered appropriate for them to inject power.

These curves describe a potential behavior and the three parameters that characterize them can be adjusted from a number of real samples or from the curve itself to obtain the continuous version. In other words, they are power functions that introduce nonlinearities to the problem and must therefore be linearized. This approximation error is not remarkable since the performance reaches high values at relatively low powers. In this case, a discretization of the curve into four intervals is chosen. In order to model the active power output, parameters are created that represent the maximum power limit that can be generated at any given time by each system on each bus and in each scenario. These parameters are characterized by the super index max. Then, for example, the maximum PV active power produced at each time *t* is given by the following formula:

$$P\_{pv}^{\text{max}}(t) = \begin{cases} \eta\_1 P\_{pv}(t) & P\_{pv}(t) \in [D\_b, P\_1] \\ \eta\_2 P\_{pv}(t) & P\_{pv}(t) \in (P\_1, P\_2] \\ \eta\_3 P\_{pv}(t) & P\_{pv}(t) \in (P\_2, P\_3] \\ \vdots \\ \eta\_t P\_{pv}(t) & P\_{pv}(t) \in (P\_{u-1}, P\_u] \end{cases} \\ = \begin{cases} 0.98 P\_{pv}(t) & 0.4 \le P\_{pv}(t) / P\_{nom} \\ 0.972 P\_{pv}(t) & 0.2 \le P\_{pv}(t) / P\_{nom} < 0.4 \\ 0.955 P\_{pv}(t) & 0.08 \le P\_{pv}(t) / P\_{nom} < 0.2 \\ 0 & \text{any other case} \end{cases} \tag{40}$$

where *Ppv*(*t*) is the active power coming from the solar field at time *t*, *η<sup>n</sup>* are the yields corresponding to the average of the discretized interval, *Pnom* is the nominal power of the inverter of the solar field and *Pu*, *Pu*−<sup>1</sup> are the powers of each interval. This expression is applied to each PV system in each bus and scenario. In the specific case of the wind turbine with doubly fed induction machine (DFIG) this efficiency only applies to the percentage of the power that circulates through the converter, that is, about 30% of the generated power.

Another important component of these machines is the ability to absorb/inject reactive power. Based on the specific characteristics of each manufacturer's model and for each application, the operating region of the converters is limited differently. In the case of PV inverters, the active power is limited to non-negative values so that the portion of the circle with the positive semi-axis of abscissa (right half of the circle) is obtained, in other words, the following constraints apply:

$$\begin{bmatrix} \left[ \sin \left( \frac{2\pi l}{k} \right) - \sin \left( \frac{2\pi (l-1)}{k} \right) \right] p\_{\text{pv}}(t) - \left[ \cos \left( \frac{2\pi l}{k} \right) - \cos \left( \frac{2\pi (l-1)}{k} \right) \right] q\_{\text{pv}}(t) \le s\_{\text{inv}} \sin \left( \frac{2\pi}{k} \right) \\\\ 0 \le p\_{\text{pv}}(t) \le P\_{\text{pv}}^{\text{max}}(t) / S\_{\text{base}}, \quad -s\_{\text{inv}} \le q\_{\text{pv}}(t) \le s\_{\text{inv}} \end{bmatrix} \tag{41}$$

where *ppv* is the decision variable of active power per unit that the inverter or group of PV inverters deliver to the grid, *qpv* is the variable of reactive power per unit that the PV inverter delivers to the grid and *sinv* is the parameter corresponding to the nominal or maximum apparent power of the inverter also expressed in per unit (pu). Here *sinv* is a deterministic parameter while *ppv* and *qpv* are stochastic variables. Again, the above expression is applied to each bus where there is a PV installation.

For wind turbines, the constraint is similar except that the power factor *fpwt*, both inductive and capacitive, is limited to maximum values of 0.95, as indicated by the manufacturer [41]:

$$\begin{cases} \left[ \sin\left(\frac{2\pi l}{k}\right) - \sin\left(\frac{2\pi (l-1)}{k}\right) \right] p\_{\text{wt}}(t) - \left[ \cos\left(\frac{2\pi l}{k}\right) - \cos\left(\frac{2\pi (l-1)}{k}\right) \right] q\_{\text{wt}}(t) \le s\_{\text{wt}} \sin\left(\frac{2\pi}{k}\right) \\\\ 0 \le p\_{\text{wt}}(t) \le P\_{\text{wt}}(t) / S\_{\text{base}} - \tan(\cos^{-1}(fp\_{\text{wt}})) s\_{\text{wt}} \le q\_{\text{wt}}(t) \le \tan(\cos^{-1}(fp\_{\text{wt}})) s\_{\text{wt}} \end{cases} \tag{42}$$

The following consideration should be noted: inverters only provide/absorb reactive power from the grid if there is active power available in the solar field. Likewise, the DFIG only manages reactive power when it is possible to generate active power. In principle, this limitation is not due to the design of the inverters, but is normal practice in this type of installation. In algebraic terms, this is achieved by modifying the restriction of the reactive power limits. The photovoltaic case is shown as an example, but for all other converters it is analogous:

$$\begin{cases} -s\_{\rm inv} & P\_{pv}^{\max}(t) > 0 \\ 0 & P\_{pv}^{\max}(t) = 0 \end{cases} \le q\_{pv}(t) \le \begin{cases} s\_{\rm inv} & P\_{pv}^{\max}(t) > 0 \\ 0 & P\_{pv}^{\max}(t) = 0 \end{cases} \tag{43}$$

The converters used in storage systems are a special case since the active power flow can be both positive and negative (discharge and load), however, the approach applied is to separate these two processes (to save binary variables and thus speed up the resolution) and apply a similar constraint to each one together with the additional reactive limitations presented by these converters:

$$
\begin{bmatrix}
\sin\left(\frac{2\pi l}{k}\right) - \sin\left(\frac{2\pi (l-1)}{k}\right) \\
\sin\left(\frac{2\pi l}{k}\right) - \sin\left(\frac{2\pi (l-1)}{k}\right) \\
\end{bmatrix}
\begin{bmatrix}
p\_{ESS2\text{net}}(t) - \left[\cos\left(\frac{2\pi l}{k}\right) - \cos\left(\frac{2\pi (l-1)}{k}\right)\right] q\_{ESS}(t) \le s\_{ESS} \sin\left(\frac{2\pi}{k}\right) \\
\cos\left(\frac{2\pi l}{k}\right) - \cos\left(\frac{2\pi l}{k}\right) - \cos\left(\frac{2\pi (l-1)}{k}\right) \left[q\_{ESS}(t) \le s\_{ESS} \sin\left(\frac{2\pi}{k}\right)
\end{bmatrix}
\tag{44}
$$

where the subscript *ESS*2*net* indicates power transfer from the storage system to the grid and *net*2*ESS* from the grid to the storage system. These and the following constraints are for each *ESS,* at each bus, at each time and in each scenario. The capacitive (*cap*) and inductive (*ind*) reactive limits are defined by giving a range to the variable as follows:

$$-\sin(\cos^{-1}(fp\_{\text{ESS},\text{ind}}))s\_{\text{ESS}} \le q\_{\text{ESS}}(t) \le \sin(\cos^{-1}(fp\_{\text{ESS},\text{cap}}))s\_{\text{ESS}}\tag{45}$$

where *sESS* is the nominal apparent power of the converter of each storage system.

In this way, the P-Q operating curve of these converters is approximated in such a way that there is a circle cut by a horizontal straight line at the top (capacitive power factor) and by another horizontal straight line at the bottom (inductive power factor). This approximation is conservative in the sense that it does not consider power peaks slightly above nominal, as these converters allow according to the manufacturer, nor does it represent the performance variation in a non-linear way in two quadrants: first and fourth. Note that these considerations could be added by combining circles that cut in such a way that the operating area is only that enclosed by all of them at the same time, however, this increases the computational time and it is not considered meaningful to implement this detail in this case.

The complete energy balance corresponding to the pu power generation can now be expressed as shown below:

$$p\_n^G(t) = p\_{wt}^n(t) + p\_{pv}^n(t) + p\_{ESS2net}^n(t) - p\_{net2ESS}^n(t) \tag{46}$$

$$q\_n^G(t) = q\_{wt}^n(t) + q\_{pv}^n(t) + q\_{ESS2net}^n(t) - q\_{net2ESS}^n(t) \tag{47}$$

Here all variables are expressed for each time point *t*, each bus *n* and stochastic scenario.

#### *4.5. Energy Consumption*

The consumption data are stochastic parameters as they vary over time, for each bus and each scenario. They are defined on the basis of the demand curve proposed in the literature [42,43], but modified to make it similar in shape to the Spanish typical curve. In turn, in the aforementioned literature, it is proposed that different loads depend on whether they correspond to the residential or industrial sector, and we have sought to maintain this distinction. The power factor of each bus is also known from the benchmark and the nominal power is given in its apparent form. The network demand is given as a percentage of the nominal power of each bus. Then, the active power to be supplied at each bus n is given by the following expression:

$$p\_n^D(t) = \frac{1}{S\_{base}} (S\_{indus}(n) f p\_{indus}(n) rate\_{indus}(t) + S\_{res}(n) f p\_{res}(n) rate\_{res}(t)) \tag{48}$$

where *res* denotes residential and *indus* denotes industrial, furthermore here *S* corresponds to the nominal apparent power of each bus. Equivalently, the reactive demand is constructed as follows:

$$q\_n^D(t) = \frac{1}{S\_{\text{base}}} \left( \mathcal{S}\_{\text{indus}}(n) \sqrt{1 - f p\_{\text{indus}}(n)^2} rate\_{\text{indus}}(t) + \mathcal{S}\_{\text{res}}(n) \sqrt{1 - f p\_{\text{res}}(n)^2} rate\_{\text{res}}(t) \right) \tag{49}$$

#### *4.6. CO2 Emissions and Other Factors*

Finally, to define the objective function, it is necessary to estimate the CO2 emissions corresponding to the import of energy from the transmission grid. This is achieved by applying a time-varying emissions factor that represents the number of metric tons of CO2 that it costs to produce one MWh unit of energy. This factor is considered to be similar to the Spanish factor and depending on the case study will be constant over time or may vary. This factor is deterministic because there is no correlation with other variables and there are no studies of its prediction.

The CO2 emissions corresponding to the import of energy from the slack bus over a period T are calculated using the following expression:

$$Emission = \sum\_{t}^{T} f\_{\text{CO}\_2}(t) p\_{\text{slack}}(t) S\_{\text{base}} \Delta t \tag{50}$$

The effect of this factor is studied in different sections. Note that the slack bus powers are free variables and if, for example, energy export is not allowed or the power factor is to be limited, constraints such as the following must be added:

$$p\_{slack}(t) \ge 0, \ q\_{slack}(t) \ge 0 \tag{51}$$

#### **5. Heuristic Methods**

Heuristic methods have been commonly used in the field of optimization, being able to obtain solutions quickly to complex problems that, even with the use of clusters of computers, may not have an optimal solution. To do so, they sacrifice solution accuracy at the cost of reduced computational cost and time. Nevertheless, they are capable of finding an exact solution to problems of relative simplicity.

These methods are often used in off-grids to find a suitable location for the development of this type of microgrids [44,45]. In recent years, research is being carried out on certain metaheuristic methods, such as Particle Swarm Optimization and Salp Swarm Algorithm, among others, which are population-based methods with promising results [46].

Heuristic methods are also applied, with success, to improve the effectiveness of the distribution network by means of reconfiguration. Reconfiguration of the distribution network aims to find the optimal combination of all switches in the distribution, mainly determining proper sizing and siting of DG together with network reconfiguration. In this way, some researchers, like Muhammad et al. [47], use the discrete network reconfiguration of the data set method, employing the Water Cycle Algorithm (WCA) together with dataset approach to reduce the complexity of search space. This type of method has good convergence performance, and can obtain a global optimal solution for singleobjective optimization problems. Others, such as Helmi et al. [48], propose the Harris Hawks Optimization (HHO) to minimize the power losses of the network.

However, all heuristic methods that are inspired by natural processes have parameters that are highly dependent on their own algorithm; therefore, the algorithm may behave differently affecting its performance [49]. Furthermore, according to Yang, and in relation to computational cost, no consensus has been reached on what are the best values or configurations of an algorithm, nor on possible ways to adjust these parameters to achieve the best performance [50]. On the other hand, the selection of a method as the most appropriate for solving a problem such as energy management in hybrid microgrids is an open problem [51]. This is largely based on the "no free lunch" theorem [52] of mathematical optimization, which shows that at the same time that a heuristic is very efficient for one collection of problems, it is very inefficient for another collection.

So, there are a multitude of methods that can be applied in this field, without any of them being clearly better than any other for any topology and type of grid, as we have discussed in previous sections. Even though the optimal power flow is a non-convex problem [53], convex approximations for the power flow equations have been studied [54,55], but generally assuming strong approximations, such as all generators are constant current injections, which is far from real microgrids and further from hybrid microgrids [56].

In this paper, we try to optimize energy management without reconfiguring the network, assuming that the location and schedule of DG and storage banks does not depend on the utility, but on the consumers, and therefore, it is not in their hands to reconfigure the network. Thus, we will not take into account, for example, the temporary shifting of loads, which is currently possible with electric vehicles. The study of controllable loads could be attacked by stochastic optimization algorithms, according to Hosseini et al. [57] or Barbato et al. [58]. In the following, two heuristic methods known in the study of microgrids will be described, namely: Hysteresis Band Control and Fuzzy Logic Control.

#### *5.1. Hysteresis Band Control*

The heuristic method of Hysteresis Band Control has been a certain success in the field of microgrids. According to the example proposed by Ipsakis et al. [59], from which

this control strategy will be explained, a hysteresis band consists of the operating range existing between two limit values of a problem variable, in this case, the State of Charge (SOC) or State of Charge of the accumulator. In this way, the storage units of the microgrid in the example, electrolyzer and fuel cell, absorb energy or give it up according to the band defined by the limit values mentioned:


This way of operating has benefits such as reducing the number of start-ups and shutdowns of the electrolyzer, which can reduce its life expectancy. Also, the example method achieves a protection of the accumulator by excluding it from excessively long operation.

Appropriate SOC limit values must be calculated to achieve effective control and in accordance with the operating requirements of the problem.

#### *5.2. Control by Means of Fuzzy Logic*

The Fuzzy Logic Control (FLC) method is characterized by reaching solutions to problems in which the data, variables, or in short, the available information, is ambiguous or imprecise, hence the term fuzzy. Therefore, fuzzy control is characterized by being described by what could be called discard logic, that is: If a certain event occurs, then the control signal will take the value 'X'.

This is why Fuzzy Logic Control (FLC) lacks accuracy when it comes to providing robust solutions. However, Fuzzy Logic Control has some advantages that make it attractive for tackling certain types of problems, these are [60,61]:


We can brief the concept by explaining its methodology with the following steps: firstly, the input data provided are processed and a smearing or merging is performed on them. In this first step, certain qualitative characteristics are given a numerical value. Once this is done, decisions are made in accordance with logical relationships called Fuzzy Rules. Finally, the defuzzification process takes place, in which concrete data are obtained that will be used to generate the appropriate control signals required by the problem.

Despite the ease of implementation of the heuristic methods described above, the large number of restrictions and variables that appear in the microgrid under study make them a bad strategy to follow for its control, as it is difficult to find optimal solutions to the problem [62,63].

#### **6. Machine Learning Methods**

Machine Learning is a scientific discipline that tries to make systems learn automatically. Learning, in this context, means identifying complex patterns in millions of pieces of data [64,65]. The machine that actually learns is an algorithm that reviews the data and is able to predict future behavior in some fields of knowledge [66–68]. Machine Learning is therefore a process of knowledge induction, that is, a method of deriving a general statement by generalizing from statements describing particular cases.

Machine Learning is learning from data, it is discovering the structure and patterns underlying the data. The main objective of Machine Learning is to extract the information contained in a dataset to acquire knowledge to make decisions about new datasets [69]. Formally, and according to Mitchell [70], we can define the algorithms used by Machine Learning as:

"A computer program is said to learn from experience *E* with respect to some class of tasks *T* and performance measure *P* if its performance at tasks in *T*, as measured by *P*, improves with experience *E*".

These learning algorithms are based on a set of data on which to learn and then apply the experience gained on other sets. It is necessary to evaluate its performance on a set other than the one on which the system has been trained in order to obtain a valid estimate of its generalizability to new examples. Thus, the available data set is divided into two subsets: on the one hand, we have the training set and, on the other one, the validation set or test set. In this way, the model is generated from the training data and evaluated in the test set, in which the accuracy of the model can be measured. The obtained result on this set is a good approximation to the expected one for the new data [71].

Therefore, generalization is one of the key aspects in the design of Machine Learning algorithms [72]. At the same time, the models must fit the training set and capture all its information. In this way, the problem of balancing bias and variance arises: bias measures the average error of the model using different training sets, while variance measures the sensitivity of the model to small changes in the training data [73]. In other words, very complex models have a low bias and a high variance, which is known as overfitting. On the other hand, simple models have a high bias but a very low variance. Overfitting occurs when, by adding levels to the Decision Tree, the hypotheses are so refined that they describe the examples used in the learning process very well; however, when evaluating the examples, an error occurs. That is, it classifies the training data very well, but then it fails to generalize the test set. This is because it learns down to the noise of the training set, adapting to the regularities of the training set [74]. Therefore, overfitting will be an important evaluation indicator to take into account in the study.

Machine Learning algorithms are usually divided into three categories, the first two being the most common:


#### *6.1. Operation with Machine Learning Models*

The process to be followed for the construction of a Machine Learning system can be divided into:

1. Data collection. This is usually a tedious process that takes up a large part of the development of the system, since it is generally necessary to collect large amounts of data in order to ensure that the used sample is representative of the set under study.

2. Feature selection. This is a critical step since it is necessary to extract those variables that are useful to distinguish the patterns of each category.

3. Choice of model. In this step, we will choose the model that best fits our problem and that achieves the expected performance on the test set. This model, among other tasks, must maintain the bias-variance balance explained above.

4. Model training. In this phase, the classifier is built, whose parameters are adjusted from the training data set. Finding the parameters that fit our model is an optimization problem since the objective is always to minimize a certain objective function.

5. Evaluation of the model. Using the test set, an error measure is set and the performance of the model is obtained. If the result is not expected, it is necessary to test by going back to each of the previous points and go through the process again.

In mathematical terms, the principle of Machine Learning, in a supervised learning context, consists of starting from a sample of learning:

$$\mathcal{L} = \left\{ (\mathbf{x}\_{\text{tr}} y\_{\text{n}}) | n = 1, 2, \dots, N, \mathbf{x}\_{\text{n}} \in \mathbb{R}^{d}, y\_{\text{n}} \in \{1, 2, \dots, \mathcal{C}\} \right\} \tag{52}$$

constituted by *n* realizations of a pair of random variables (*X*, *Y*), to construct a *<sup>f</sup>* : *<sup>d</sup>* <sup>→</sup> {1, 2, . . . , *<sup>C</sup>*} function which, given a new *<sup>X</sup>* input vector, can predict with some degree of certainty the variable *Y* = *f*(*X*). For each observation (*xi*, *yi*) of L, the variable *xi* ∈ *X* is called the input variable or explanatory variable and *yi* ∈ *Y* the dependent or output variable [75]. When the dependent variable is discrete or categorical, we speak about a classification problem; and when it is continuous, about a regression problem. That is to say, depending on the type of objects that we are trying to predict, there are two types of problems:


One of the most widely used Machine Learning methods is Decision Tree Learning. This is a method for approximation of discrete-valued functions, robust to noisy data and able to learn disjoint expressions. There is a family of Decision Tree Learning algorithms: ID3, C4.5, ... In turn, based on these Decision Trees, hybrid methods have been created that build more than one Decision Tree: Bagging, Boosting and Random Forest.

#### *6.2. Decision Trees*

Learning through Decision Trees is based on the principle of divide and conquer. Let L a sample be defined as:

$$\mathcal{L} = \left\{ (\mathbf{x}\_{n\prime} y\_n) | n = 1, 2, \dots, N, \mathbf{x}\_n \in \mathbb{R}^d, y\_n \in \{1, 2, \dots, \mathbb{C}\} \right\} \tag{53}$$

where *N* is the number of elements in the data set, *C* the number of distinct classes and *d* the number of variables defining the examples *xn* in the set. Each of these examples is represented by a vector *xn*, which has its corresponding class label *yn* associated with it and it is defined by different variables, which can be numerical (their values are real numbers) or categorical (they take values in a finite set in which there is no ordering relationship). Sometimes, these vectors *xn* are also cited as feature vectors.

A *T* Decision Tree is an ordered sequence of decisions, normally from questions, in which the next decision depends on the answer to the current one. These decisions are taken normally from questions that are formulated on the variables that define each *x* element in order to assign them a *y* certain class. This process, including its corresponding questions, decisions and bifurcations, is naturally represented by means of a tree.

In a Decision Tree, each node of the tree is an attribute (field) of the examples, and each branch represents a possible value of that attribute. The first node is known as the root node, which is successively connected to the other nodes until it reaches the leaf nodes, those that have no descendants, that is, the end of the branches of the tree. Each node is assigned one of the questions of the sequence, while each leaf node is assigned a class label.

In this way, the question of the root node is asked of the whole set L, which is subdivided until reaching the last nodes (leaves), which constitute a disjoint partition of the initial feature space. This happens because, when given a node, one and only one branch will be followed by each instance of the training set.

Decision Trees perform well with large volumes of data, as it does not require loading all data into memory at once. The computation time scales well with a linearly increasing number of columns [76].

The advantages of Decision Trees are that they are easy to understand and interpret, rule generation is simple, it reduces the complexity of the problem, and training time is not very long.

One of the disadvantages is that if an error is made at a high level, successive nodes would be poorly created. In the construction of a Decision Tree, the most complicated step is to determine which attribute to base a node on, because if there are many features, the algorithm would have many options for the training data, and it would be difficult to construct a Decision Tree.

However, Decision Trees can give good results if they are combined with ensemble methods. With these methods, instead of learning a single model, several models are learned, and the estimates from each model are combined.

Ensemble methods are combinations of models. In these techniques, it is necessary both to define how different models are to be created and how the results of each model will be combined to generate the final prediction. The aim of ensemble methods is to produce a better prediction than individual models (individual members of the ensemble).

The most common ensemble methods are Bagging, Boosting and Random Forest. In all of these methods, the training set is manipulated, but in each case with a different strategy [77].

In Bagging, different samples are extracted from the training set (bootstrap samples), and these bootstrap samples are used as if they were the true training set. Boosting, on the other hand, always works with the full data set, that is, the complete dataset is always used. In Boosting, we can manipulate the weights of the data in the training set to generate different models. At each iteration, Boosting learns a model that minimizes the sum of the weights of the misclassified data.

Finally, Breiman [78] presented an ensemble method called Random Forest where bagging is used together with a random selection of attributes. At each node of each tree in the forest, a subset of the available attributes at that node is randomly selected and the best of them is selected according to the splitting criteria used in the base algorithm. The number of randomly selected attributes is an input parameter.

#### **7. Results**

We have evaluated the following of the proposed models explained in the previous sections, in order to measure the accuracy: Hysteresis Band Control, Fuzzy Logic Control, and Decision Trees (DT). We propose the IEEE microgrid test system of 69-bus [79]. The 69-bus distribution network has a nominal voltage of 12.66 kV. Its base apparent power is 10 MVA. This system has 69 nodes and 73 branches, including tie-lines, as shown in Figure 1. The order of each branch is assumed to be that of the furthest node minus one unit, except for the tie-lines in the base scheme which run from 69 to 73. Thus, a total of 73 remote switches are installed in the network, 68 of which are sectioned and ready for possible reconfiguration. As mentioned, and described by Lan et al. [80], two wind units, two solar panels and some switches are located in the microgrid, and a battery storage unit is also installed in the hybrid microgrid. We have followed the same location:

**Figure 1.** IEEE 69-bus test system.

In this article, we only show a first approximation of the methods to verify their viability and check primary results. In addition to achieving adequate energy management using all the tested models, regardless of the calculation time or the power and dedication of the used computers, we can highlight that for loads without large hourly differences (as it would have been the case for valley and peak periods with a large difference in the power value), the system can be controlled by applying the aforementioned modeling.

From the test results, it can be found the Machine Learning model (Decision Tree) can also recognize the states accurately for distribution systems, as we show in Table 1. We used the root medium square error for the comparison.

**Table 1.** Root Medium Square Error (RMSE) (voltage of the network) with different models applied to the IEE 69-bus distribution network.


HBC = Hysteresis Band Control; FL = Fuzzy Logic; DT = Decision Tree.

As we can see, the three tested methods have enough viability, taking into account the characteristics of the network. We cannot conclude a better performance for the Decision Tree-based model, because, although its result is better than the others, the test is based in a single evaluation and more evaluations with different cases (for example, with different cases of distribution of solar and wind generation and different loads) should be done in the future.

In the case of the most modern method, Decision Trees, it should be noted that both solar and wind generation have followed generation patterns established in advance, according to the data of Lan et al. [80]. Obviously, in an analysis with a real system, these future generation data should be based on historical data and take forecasting into account. In Machine Learning methods, the reliability of data and predictions is very important, as they are the foundation of this type of modeling; therefore, it is advisable to use more than one database, which will usually reduce the error of calculation and analysis. These used patterns have been the same for any unit of the same generation type; that is, all solar units follow the same pattern, and all wind units follow the same pattern. This has been done to simplify the performed analysis and in accordance with Kovousi–Fard and Khodei [81].

The evaluations are based on the values shown in Figure 2 (base case), which represents the total energy consumption over a standard day.

**Figure 2.** Total energy consumption value (kWh) over a standard day.

Figure 3 shows the values obtained in the energy storage units, including the battery allocated at node 15. It presents the 24-h values of the state of the charge (SOC) of all the devices (that is, those related to the PV systems and the aforementioned battery unit). As can be seen, the energy storage units are strongly influenced by the solar behavior and are discharged from 16–17 h onwards, as the PV power generated gradually decreases.

**Figure 3.** 24 h SOC (%) of three of the energy storage devices of the microgrid, allocated at node 18 (orange); node 28 (green); node 56 (blue).

Figures 4 and 5 show the two main cases studied: Figure 4 corresponds to the situation where there is no distributed generation, and Figure 5 where there is distributed generation. As can be seen in Figure 4, the voltage profile has been improved with the use of the three methods, not being able to conclude which of the three is the best, mainly because, although the three improve the initial case, depending on the section of the system to be analyzed, the best method is one or the other. The tie-switch distributions generated by each method are shown in Table 2.

**Figure 4.** Case of no distributed generation. Voltage profile. Base case (blue); Hysteresis Band Control method (orange); control by means of Fuzzy Logic (grey); Decision Trees method (yellow).

**Figure 5.** Case of distributed generation. Voltage profile. Base case (blue); Hysteresis Band Control method (orange); control by means of Fuzzy Logic (grey); Decision Trees method (yellow).

**Table 2.** Case without distributed generation: tie-switches of the configuration of the system.


The IEEE-69 system has also been tested with distributed generation. As mentioned above, the distributions of wind and solar generation have followed the one published in the literature references, but Table 3 shows the averages of a typical day for the four nodes where they are located (6 and 68 for wind generation, and 25 and 50 for solar generation). In addition, Table 3 shows the tie-switches for each method. It can be seen in Figure 5 how, in general, the three methods improve the voltage profile of the system again after the reconfiguration of the system, being slightly different the solution (topology) found by each of them.


**Table 3.** Case with distributed generation: tie-switches of the configuration of the system and power (node–kW).

HBC = Hysteresis Band Control; FL = Fuzzy Logic; DT = Decision Tree.

#### **8. Conclusions**

In this article, we presented the optimization of the operation of electrical hybrid microgrids, focused particularly on the mathematical modeling. The set of loads (consume), generators (with a predominance of the renewable energies) and storage systems (at the present, particularly batteries) makes up the electrical hybrid microgrid. The management of this type of networks presents many challenges and various options to be implemented.

The mathematical modeling of the different components is a very important step in the control and management that, in recent years, is increasing with the irruption into the market of new technologies or new models for the use of energy. The main contribution is the mathematical modeling of several components of the hybrid microgrid. This modeling can be used in different methods of control and management of the network, and its feasibility has been shown in three methods, including one based on the Decision Tree method, which belongs to the Machine Learning family. The results on a test system of 69 buses show that its implementation is possible.

The three methods have been compared both in the case of existence of distributed generation and in the case of its non-existence, in order to obtain a better view of their behavior. The three methods improve the voltage profile of the system, using a different topology, although similar (the tie-switches used are topologically very close or even the same), demonstrating their effectiveness. The integration of distributed generation in the system causes small differences in the results of each method, these differences being more noticeable the longer the line (or branch) in the system.

Nevertheless, although the results indicate that the Decision Trees method is partially better than other algorithms, more tests are needed and they need to be carried out with different typologies, not only of the network itself but also of its components. As a future research line, further tests with different generation and load levels can be identified.

**Author Contributions:** Conceptualization, J.B. and O.G.; methodology, J.B., O.G., C.R. and C.V.; formal analysis, J.B.; investigation, J.B., E.B. and O.G.; writing—original draft preparation, J.B.; writing—review and editing, J.B., E.B. and C.R.; supervision, J.B.; project administration, J.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We thank the support of the laboratory of the Applied Mathematics Department of the University of the Basque Country, UPV/EHU, for its technical support.

**Conflicts of Interest:** The authors declare no conflict of interest. Authors had the role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

#### **References**

