**Advanced Operation and Maintenance in Solar Plants, Wind Farms and Microgrids**

Editors

**Luis Hern´andez-Callejo Maria del Carmen Alonso Garc´ıa Sara Gallardo Saavedra**

MDPI ' Basel ' Beijing ' Wuhan ' Barcelona ' Belgrade ' Manchester ' Tokyo ' Cluj ' Tianjin

*Editors* Luis Hernandez-Callejo ´ Ingenier´ıa Agr´ıcola y Forestal Universidad de Valladolid Valladolid Spain

Maria del Carmen Alonso Garc´ıa Energy CIEMAT Madrid Spain

Sara Gallardo Saavedra Ingenier´ıa Agr´ıcola y Forestal Universidad de Valladolid Valladolid Spain

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: www.mdpi.com/journal/applsci/special issues/ advanced renewables).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4062-7 (Hbk) ISBN 978-3-0365-4061-0 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


Reprinted from: *Appl. Sci.* **2021**, *11*, 11332, doi:10.3390/app112311332 . . . . . . . . . . . . . . . . **153**


## **About the Editors**

#### **Luis Hern ´andez-Callejo**

Prof. Dr. Luis Hernandez-Callejo is Electrical Engineer at Universidad Nacional de Educaci ´ on´ a Distancia (UNED, Spain), Computer Engineer at UNED and Ph.D. at Universidad de Valladolid (Spain). Professor and researcher at the Universidad de Valladolid. His areas of interest are: renewable energy, microgrids, photovoltaic energy, wind energy, smart cities, and artificial intelligence. Prof. Dr. Hernandez Callejo is an editor in numerous scientific journals, and is a ´ guest editor in many Special Issues. He has directed four Doctoral Theses, and at the moment, he is directing six Doctoral Theses. He is a professor in wind energy, solar energy, microgrids, and he collaborates with many universities in Spain and in the rest of the world.

#### **Maria del Carmen Alonso Garc´ıa**

Dr. Carmen Alonso-Garc´ıa is a Staff Scientist in the Photovotlaic Solar Energy Unit of Centro de Investigaciones Energeticas, Medio-ambientales y Tecnol ´ ogicas (CIEMAT) in Madrid, Spain. She ´ has performed research in photovoltaic solar energy since 1990, focusing her specialization in the characterization and modelling of PV generators, operation and degradation of different technology modules and systems, accelerated testing and circularity of the systems. She has been involved in at least 41 R+D+I projects in the field of photovoltaic solar energy, being the scientist in charge in many of them. With a number of publications in indexed journals, conferences and scientific dissemination media, her contributions cover a broad spectrum within solar photovoltaics, from basic research to applied research and demonstration projects. She has been actively involved in standardization committees and other relevant photovoltaic solar energy international groups.

#### **Sara Gallardo Saavedra**

Dr. Sara Gallardo Saavedra is a professor and researcher at the Campus Duques de Soria of the University of Valladolid (Spain). Her research focuses on the detection, characterization and classification of defects in photovoltaic (PV) modules through the use of thermography, electroluminescence, I-V curves and visual analysis. She has participated in numerous national and international R+D+I projects, carrying out an active dissemination of the results, with high regularity in the scientific production, including scientific publications in high impact factor journals, book chapters, and contributing to congresses on advanced maintenance in PV. She has made a predoctoral and a postdoctoral stay in the Unit of Solar PV Energy in the Energy Department of the Energy Research Center, Environment and Technology (CIEMAT) in Madrid and the researcher has collaborated with different institutions as University of Gavle in Sweden, with the Universidad del Valle in Colombia, with the National Polytechnic Institute of Mexico and with the University of Cuenca in Ecuador.

## **Preface to "Advanced Operation and Maintenance in Solar Plants, Wind Farms and Microgrids"**

The development of renewable generation plants and microgrids is a reality. Every day more facilities of this type are springing up, and their advance requires new research studies. Among renewable energy, solar plants (photovoltaic, thermal, and hybrid) and wind plants have had the greatest impact in recent years. With regard to microgrids, these scenarios are integrators of local generation sources (renewable or nonrenewable) and are facilities that promote energy sustainability.

In any case, both renewable generation plants and microgrids require technological development tools, mainly with regard to their operation and maintenance. Therefore, this Special Issue includes reviews, research articles, case studies, and technical notes on "Advanced Operation and Maintenance in Solar Plants, Wind Farms and Microgrids"

For solar plants, wind farms, and microgrids, the articles focus on one of the following topics:


## **Luis Hern ´andez-Callejo, Maria del Carmen Alonso Garc´ıa, and Sara Gallardo Saavedra** *Editors*

**Miguel Aybar-Mejía 1, \* , Lesyani León-Viltre 2,3 , Félix Santos 3,4 , Francisco Neves 5, \* , Víctor Alonso Gómez <sup>6</sup> and Deyslen Mariano-Hernández 1**


**Abstract:** A smart microgrid is a bidirectional electricity generation system—a type of system that is becoming more prevalent in energy production at the distribution level. Usually, these systems have intermittent renewable energy sources, e.g., solar and wind energy. These low voltage networks contribute to decongestion through the efficient use of resources within the microgrid. In this investigation, an energy management strategy and a control scheme for DG units are proposed for DC/AC microgrids. The objective is to implement these strategies in an experimental microgrid that will be developed on the INTEC university campus. After presenting the microgrid topology, the modeling and control of each subsystem and their respective converters are described. All possible operation scenarios, such as islanded or interconnected microgrids, different generation-load possibilities, and state-of-charge conditions of the battery, are verified, and a seamless transition between different operation modes is ensured. The simulation results in Matlab Simulink show how the proposed control system allows transitions between the different scenarios without severe transients in the power transfer between the microgrid and the low voltage network elements.

**Keywords:** microgrid; control system; storage system; wind turbine; primary control

#### **1. Introduction**

The constant growth in renewable energy generation and the integration of these systems into large-scale grids represents a challenge for the proper functioning of the electrical system [1]. Grid integration requirements have therefore become a significant concern as renewable energy sources, such as wind and solar photovoltaic (PV) systems, slowly begin to replace conventional plants [2]. The proper control and operation of electric microgrids integrated into the electrical network can improve the electrical system's power quality, stability, and reliability [3].

In [4], a hybrid AC/DC microgrid has been incorporated into a low voltage AC distribution system. A power control scheme is presented to improve the system's stability after load steps when the microgrid operates in island mode. However, the authors of this study do not consider the presence of mini wind power sources in the microgrid. Further, they do not describe the control algorithm for islanded operation, nor the batteries being

**Citation:** Aybar-Mejía, M.; León-Viltre, L.; Santos, F.; Neves, F.; Gómez, V.A.; Mariano-Hernández, D. Modeling and Control of a Microgrid Connected to the INTEC University Campus. *Appl. Sci.* **2021**, *11*, 11355. https://doi.org/10.3390/ app112311355

Academic Editor: Amjad Anvari-Moghaddam

Received: 4 November 2021 Accepted: 25 November 2021 Published: 30 November 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/).

fully charged—a condition that requires that distributed generators (DG) units abandon maximum power point tracking (MPPT) and share the total load in proportion to their available primary powers.

A power management system for hybrid AC/DC microgrids, designed to optimize a cost function that considers the maximum utilization of renewable resources, minimal usage of fuel-based generators, extending the lifetime of batteries, and limited utilization of the main power converter between the AC and DC micro-grids is presented in [5]. However, the power management algorithm needs several input variables which are usually not available in most microgrids, e.g., PV and wind primary energies, namely, demanded power (sum of load power, power line loss, and power loss due to the circulating current), the state-of-charge (SOC) of the battery banks, and some statistical and dynamical operational limits.

The management and optimization of a hybrid AC/DC microgrid are still open issues, as affirmed in [6]. In this study, the authors propose a model predictive power and voltage control (MPPVC) method for providing microgrid control under different possible scenarios with smooth transients. However, the precise details about the control of the DG units for off-MPPT operation while ensuring proper total load sharing are not given.

In this paper, a power management scheme for a hybrid AC/DC wind–PV–ESS microgrid is proposed. The only input variables necessary for the converters' control are those usually measured, the state of the microgrid (interconnected or not), and the battery SOC. The control schemes of the converters for all operation scenarios are described in detail. The results obtained in this paper serve as a basis for a hybrid AC/DC microgrid comprising a PV system, an energy storage system (ESS) using lithium-ion batteries, DC loads, a wind microturbine with a permanent-magnet synchronous generator system (PMSG), and AC loads that are under development as part of a project financed by the Ministerio de Educación Superior, Ciencia y Tecnología (MESCyT) of the Dominican Republic which will be installed on the campus of the Instituto Tecnológico de Santo Domingo Instituto de Santo Domingo (INTEC). The objective of the project is to implement and verify the behavior of the hybrid microgrid through the intelligent management of the distributed generation (DG) units and loads and the collection of information. While the microgrid is under construction as part of the energy management and control system (EMCS) design, its behavior is studied through a detailed simulation system. In the present study, the sizes of the photovoltaic, wind, and storage systems are the same as those of the equipment available for installation, and energy production is estimated based on forecasting data. However, implementing the control and management of this system is not limited to the specific elements of the microgrid. The results can be extrapolated to other low voltage infrastructures to meet the voltage and electrical frequency conditions. This project aims to promote the applications of hybrid DG systems integrated into electrical microgrids.

In the present work, we evaluate a control strategy that allows power transitions between the microgrid connected to the university campus and all possible operating scenarios, such as islanded or interconnected microgrids, different generation-load possibilities, and state-of-charge conditions of the battery. The converters have a low number of switches but allow adequate control of each element of the microgrid, so that the GD units operate in maximum power point tracking (MPPT) mode or promote the balance between the generated power, loads, and battery charge/discharge powers.

The overall management strategy, together with the proposal of control schemes for the DG units for off-MPPT operation while ensuring proper total load sharing, are the main original contributions of this article.

The rest of the paper is organized as follows. The proposed topology and modes of operation are described in Section 2. Models were constructed to enable the design of the system control loops, the PV, PMSG-wind, and ESS systems, together with the respective interface converters, and these are presented in Section 3. The control strategies for the interface converters are explained in Section 4. The operational scenarios used to verify and validate the proposed control strategies are described in Section 5, and the corresponding simulation results are presented in Section 6. The conclusions of the study are drawn in Section 7.

#### **2. Microgrid Topology and Modes of Operation**

Figure 1 presents a schematic diagram of the PV–wind–ESS hybrid AC/DC microgrid. The photovoltaic panels, lithium-ion battery, and PMSG-based wind system are connected to the microgrid DC bus via DC/DC converters. A bidirectional inverter is used to control the energy flow between the DC and AC microgrid buses. AC loads are connected to the microgrid at 220 V, 60 Hz AC bus, through which the microgrid is connected to the central power system. The capacity of the microgrid simulated for a PV system is 500 Wp, wind turbine capacity is 2 kW, and lithium storage capacity is 5 kWh. A simple boost converter has been chosen as the PV system interface, and a low-cost two-switch bidirectional converter topology was used for the proper control of the ESS.

**Figure 1.** Topology of the PV–wind–ESS hybrid microgrid.

The DG and ESS interface converters must allow different control modes, depending on whether the microgrid is connected to the central power system or not. Furthermore, with respect to island mode, the converters' control strategies depend on the batteries' state of charge (SoC) and the total DG units' instantaneous injected power versus total load demand. Thus, some quantities necessary for the decision about the microgrid mode of operation must be sent to a central control unit, which decides the DG unit's operation mode (under MPPT or with reduced generation) and ESS operation (charging or supplying the instantaneous power deficit).

For long-term operation, climatology information obtained for the project's location in the university was used (18◦29′15.9′′ N 69◦57′48.5′′ W).

When the microgrid is connected to the primary grid, the PMSG wind and PV systems must be controlled to maximize the use of available primary energy, i.e., MPPT is adopted. Furthermore, if the battery is not fully charged, its converter controller should determine the instantaneous power to be absorbed. Therefore, for operation in under-connected mode, the inverter between the DC and AC buses must transfer to (or absorb from) the AC microgrid bus the difference between the power generated and consumed by the battery and DC loads.

On the other hand, loads should preferably be fed by the DG sources when the microgrid is islanded. In case of insufficient generation, the batteries should complement

the power to meet the load's demand. However, if the generation exceeds the power required by the loads there are two possible situations that might result: (i) the battery is fully charged, in which case the DG systems can no longer operate at the maximum power point MPP and must share the total power consumed by the loads proportionally to their nominal powers; or (ii) the battery can be charged, allowing the DG systems to operate at the MPP, and the excess of generated power is used for battery charging.

It is essential to mention that the topology and control strategies must ensure proper operation under the different conditions described above and provide smooth transitions between operating modes. The energy management strategy described is represented in the flow diagram of Figure 2.

**Figure 2.** Microgrid energy management strategy.

#### **3. Modeling of the Microgrid**

In this section, the models of the main microgrid elements—the PMSG wind system, the PV module, and the battery—are briefly described.

#### *3.1. PV System*

The PV module can be characterized by a nonlinear I–V curve that changes according to the local temperature and irradiance conditions. Among the different proposals to accurately represent the PV module, the most used is the single-diode model, shown in Figure 3, due to its good compromise between simplicity and accuracy. Furthermore, it is possible to estimate the model's parameters based on the data provided in the panel manufacturer datasheet [7]. Table 1 describes the characteristics of the one-diode model elements.

**Figure 3.** Equivalent circuit of the photovoltaic cell.


**Table 1.** Model of photovoltaic modules.

Equations (1)–(4) demonstrate the equation of the PV panel for the light-generated photocurrent as presented in [7]. The output current of the single diode model is a function of the output voltage, and it can be described as:

$$I = I\_{\mathcal{S}} - I\_{\text{sat}} \left[ e^{\left(\frac{V + IR\_{\mathcal{S}}}{Vt}\right)} - 1 \right] - \frac{V + IR\_{\mathcal{S}}}{R\_p} \tag{1}$$

$$V\_t = \frac{N\_s Akt}{q} \dots \tag{2}$$

where *V* and *I* are the module output voltage and current, *I<sup>g</sup>* is the photogenerated current, *Isat* is the reverse saturation current of the diode, *V<sup>t</sup>* is the thermal voltage, *q* is the electron charge, *A* is the ideality factor of the diode, *k* is the Boltzmann constant, *T* is the module temperature, *N<sup>s</sup>* is the number of series-connected cells forming the PV module, and *R<sup>s</sup>* and *R<sup>p</sup>* are the series and the parallel resistances, respectively.

The PV module manufacturers do not directly provide the five parameters of the single-diode model. For this reason, it is difficult to determine those parameters using simple analytical methods. Instead, all datasheets provide the following information for the standard test conditions (STC): open-circuit voltage (Voc), short-circuit current (Isc), MPP voltage (Vmp), MPP current (Imp), a temperature coefficient for Voc (kV), a temperature coefficient for Isc (ki), and maximum power (Pmp). An analysis of the electrical model in three operation points of the I–V curve (Isc, Voc, and Pmp) allows the unknown parameters of the electrical model to be related to the datasheet information.

Several authors propose novel photovoltaic parameter estimation methods. This article uses the method proposed in [7], which involves finding the five unknown parameters that guarantee the absolute minimum error between the P–V curves generated by the electrical model and the P–V curves provided by the manufacturers' datasheets for different external conditions, such as temperature and irradiance. *I<sup>g</sup>* is given by (3), *Isat* is obtained from (4):

$$I\_{\mathcal{S}} = \left[I\_{\rm sc,STC} + k\_i(T - T\_R)\right] \frac{\mathcal{S}}{1000}.\tag{3}$$

$$I\_{\rm sat} = \frac{I\_{\rm g} - \frac{V\_{\rm gc}}{R\_p}}{e^{\frac{V\_{\rm g}}{\Delta T}} - 1} \dots \tag{4}$$

A complete scan of all possible values of A (from 1 to 2, with a step of 0.01) and *R<sup>s</sup>* (from 0 to 2 Ω, with a step of 1 mΩ) is made.

#### *3.2. Wind Turbine and Permanent Magnet Synchronous Generator*

The instantaneous power delivered to the PMSG axis by the wind turbine, as in the system described in [8], can be represented by:

$$P = \frac{1}{2}\rho A V^3 \mathbb{C}\_p(\lambda, \beta),\tag{5}$$

where *ρ* is the air density, *A* is the area swept by the turbine rotor blades, *V* is the instantaneous wind speed and *Cp*(*λ*, *β*) is the efficiency power conversion factor, which is a nonlinear function of the tip speed ratio *λ* = *ωtR*/*V* and the pitch blades angle *β*. Since the manufacturers of wind turbines do not give information about power conversion factor characteristics, some typical empirical curves are often used for representing wind turbines in power system studies, such as the one presented in [8]:

$$\mathcal{C}\_p(\lambda, \beta) = 0.22 \left( \frac{116}{\lambda\_i} - 0.4\beta - 5 \right) e^{\frac{-12.5}{\lambda\_i}},\tag{6}$$

where:

$$\frac{1}{\lambda\_i} = \frac{1}{\lambda + 0.08\beta} - \frac{0.035}{\beta^3 + 1}. \tag{7}$$

The PMSG was represented using the model in the *dq* reference frame rotating at a rotor electrical angular speed *ω<sup>r</sup>* = (*P*/2)*ω<sup>t</sup>* :

$$\begin{aligned} v\_{sd} &= R\_s i\_{sd} + \frac{d}{dt} \lambda\_{sd} - \omega\_r \lambda\_{sq} \\\\ v\_{sq} &= R\_s i\_{sq} + \frac{d}{dt} \lambda\_{sq} + \omega\_r \lambda\_{sd} \\\\ \frac{2I}{P} \frac{d}{dt} \omega\_r &= T\_{\ell} - T\_m - \frac{2b}{P} \omega\_r \end{aligned} \tag{8}$$

where *vsd*, *vsq*, *isd*, *isq*, *λsd*, and *λsq* are the direct and quadrature axes components of the stator voltage, current, and flux space vectors, respectively. *T<sup>e</sup>* and *T<sup>m</sup>* are the electromagnetic and mechanical torque machines, and *J* and *b* are rotor inertia and viscous friction coefficients. The flux–current relations and the electromagnetic torque equation complete the model:

$$\lambda\_{\rm sd} = L\_{\rm sd} i\_{\rm sd} + \Lambda$$

$$\lambda\_{\rm sq} = L\_{\rm sq} i\_{\rm sq} \tag{9}$$

$$T\_{\varepsilon} = \frac{3}{2} \frac{p}{2} \left(\lambda\_{\rm sd} i\_{\rm sq} - \lambda\_{\rm sq} i\_{\rm sd}\right) \\ = \frac{3}{2} \frac{p}{2} \left[\Lambda i\_{\rm sq} + \left(L\_{\rm sd} - L\_{\rm sq}\right) i\_{\rm sd} i\_{\rm sq}\right]$$
 $\text{commutment magnet flux}$ 

where Λ is the permanent magnet flux.

#### *3.3. Battery Model*

Mathematical modeling and dynamic simulation of battery storage systems can be challenging due to their nonlinear nature. The published literature has presented several thermal models of lithium-ion battery packs [9,10].

In [9], a suitable, convenient dynamic battery model is presented that can be used to model a general battery storage system. The proposed dynamic battery model can analyze the effect of temperature, cyclic charging/discharging, and voltage stabilization effects. Temperature makes a difference to three battery parameters in an equivalent electrical circuit battery model. These are the polarization voltage *K*, the battery constant *E*0, and the exponential coefficients, *A* and *B*. The dynamic battery model can be described as [9]:

$$E\_{(q,t)} = X\_{\rm E0}(T).E\_0 - X\_{\rm K}K\left(\frac{\mathcal{Q}}{\mathcal{Q}-q}\right) + A\exp^{-X\_{\rm B}B\_q} \dots \tag{10}$$

where *E*(*q*,*t*) = no load voltage (V), *E*<sup>0</sup> = battery constant voltage (V), *K* = polarization voltage (V), *Q* = battery capacity (Ah), *A* and *B* = exponential constants, and *q* = charge or extracted capacity (Ah).

The thermal effect on *A* has been ignored to make the model simpler, since *A* and *B* are highly related. For this reason, the mathematical relationships are:

$$X\_{\text{ll}} = f(T) = A + BT + CT^2 \dots \tag{11}$$

$$
\mathfrak{n} = f(\mathcal{K}, \mathcal{E}\_0, \mathcal{B}.\mathcal{C}) \tag{12}
$$

#### **4. Control Strategies**

The converters of the microgrid DG and battery units and the converter between the DC and AC buses must be controlled according to the modes of operation described in Section 2. The control scheme implemented in each converter is presented in this section. Procedures for tuning these controllers are also explained.

#### *4.1. PV Converter Control*

As already mentioned in Section 2, the PV generation system should operate in MPPT mode in three possible situations: (i) when the microgrid is in the connected mode; (ii) when the microgrid is islanded, and the total power demanded by the loads is greater than the total power produced by the DG units (the battery supplies the difference); and (iii) when the microgrid is islanded, and the power produced by the DG units is greater than the power demanded by the loads but the battery is not fully charged (so it can absorb the exceeding generated power). However, if the microgrid is islanded, the DG units' produced power is more significant than that demanded by the loads, and the battery is fully charged, then the MPPT cannot be imposed, and the generated power must be reduced so that it equals the total load power. In this last case, the load power is shared among the DG units proportionally to the DG units' rated powers.

#### 4.1.1. PV Converter MPPT Method

The electrical power supplied by PV cells is a non-linear function of voltage, current, temperature, and solar irradiance. This non-linearity makes it challenging to obtain the operating point at which its maximum power is extracted, as this point varies throughout the day due to variations in irradiance and temperature. The objective of any MPPT algorithm is to find the voltage to be applied to the PV string terminals that results in the maximum possible generated power for the momentary conditions of irradiance and temperature.

Several techniques for determining MPP have been proposed over the years. These techniques vary in complexity, speed of convergence, voltage fluctuation around the MPP, and computational cost. The MPPT technique chosen for this project was the very popular perturb and observe (P&O) algorithm due to its simplicity and acceptable performance in most situations [11]. The P&O MPPT technique changes the PV array voltage in one direction and observes the variation in the output power. If the generated power increases, then the voltage perturbation continues in the same direction. However, if the output power reduces, the voltage perturbation occurs in the opposite direction. The voltage disturbance process is repeated periodically, each *TMPPT*, and the array voltage oscillates around the MPP. This oscillation can be minimized by reducing the size of the disturbance, but minimal disturbances make the technique slow to track the MPP. The parameters of the P&O MPPT method are *TMPPT* and the magnitude of the voltage perturbation is ∆*V*. Typical choices of these parameters are ∆*V* = 0.5% of the PV array open-circuit voltage and *TMPPT* = 4*τVdc*, where *τVdc* is the time constant of the PV array output voltage variation.

After the reference value of the PV array, DC voltage is determined, and the DC/DC boost converter synthesizes it. Thus, the duty cycle of the boost converter is determined from the DC bus voltage of the microgrid and the reference PV array voltage, calculated according to the P&O algorithm for achieving MPPT.

4.1.2. PV Converter Control for Sharing the Microgrid Load Power

As has already been mentioned, if the microgrid is islanded, the battery is fully charged, and the total DG power exceeds the microgrid AC + DC load power, then each DG source must reduce the generated power until the load power demand is shared among the DG units. Ideally, this sharing should be proportional to the maximum power delivered by each DG unit.

If the DG power is greater than the load power demand, considering that no power can be transferred to the grid or the battery, the excess power flows to the DC bus capacitors bank, increasing the DC bus voltage. Thus, based on the information that the microgrid is islanded and that the battery state-of-charge is complete, the DC bus voltage increase indicates that the power generated by each DG unit must be reduced.

A typical voltage–power curve for a PV system is shown in Figure 4. However, it can be observed that, below the MPP voltage, the voltage–power relation can be considered approximately linear.

**Figure 4.** Typical P–V characteristic of a PV string.

The strategy proposed here to make the DG units share the total power required by the microgrid loads is based on the DC bus voltage increase per unit of the maximum allowed voltage increase, expressed as:

$$
\Delta V\_{\rm DC,pu} = \frac{\Delta V\_{\rm DC}}{\Delta V\_{\rm DC}^{MAX}} = \frac{V\_{\rm DC} - V\_{\rm DC}^{rated}}{V\_{\rm DC}^{MAX} - V\_{\rm DC}^{rated}}.\tag{13}
$$

Assuming that on the left side of the P–V characteristic, the PV-generated power is approximately proportional to the PV string output voltage, this voltage can be reduced linearly with the DC bus voltage increase. Since the DC bus voltage becomes constant only when the generated and load powers are equal, the steady-state condition ensures that, for some DC bus acceptable overvoltage, the PV output voltage ensures DG and load powers to match. The PV output voltage command can then be obtained, as shown in Figure 5.

If other PV units were connected to the microgrid, their output voltages would be calculated using the scheme presented in Figure 5. Thus, the generated powers would be approximately proportional to each PV unit's maximum power available.

**Figure 5.** PV output voltage references calculation for MPPT and for reduced generation to track load power demand.

#### *4.2. PMSG-Based Wind Turbine Control*

#### 4.2.1. Control for Ensuring MPPT

The PMSG stator terminals are connected to a three-leg AC/DC voltage source converter. If the PMSG rated stator voltage is low enough, the AC/DC converter can directly connect to the microgrid DC bus, otherwise a transformer between the PMSG terminals and the AC/DC converter or a DC/DC converter between the AC/DC converter and the microgrid DC bus would be necessary.

The main objective of the AC/DC converter is to optimize the power extracted from the turbine for any incoming wind speed. The wind microturbine has fixed-pitch blades at angle *β*. Thus, according to Equations (5)–(7), the turbine efficiency factor *Cp*(*λ*, *β*) depends only on the tip speed ratio *λ* = (*ωtR*)/*V*. For each wind speed *V*, there is an angular speed of the turbine rotor *ω<sup>t</sup>* that maximizes *Cp*(*λ*, *β*), allowing the maximum available power to be extracted. This optimum angular speed is obtained from:

$$
\omega\_t^{opt} = \frac{\lambda^{opt} V}{R} \tag{14}
$$

The PMSG reference rotor speed (in electrical radians per second) can then be obtained from *ω*∗ *<sup>r</sup>* = (*P*/2)*ω opt t* . A vector control scheme is then applied to impose the reference rotor speed.

The PMSG model is in a rotor-oriented reference frame, i.e., the d-axis is aligned with the permanent magnet flux. The reference d-axis stator current is then set to zero since a negative value would reduce the main flux and a positive value could cause magnetic saturation. Furthermore, if *isd* = 0, the electromagnetic torque becomes *T<sup>e</sup>* = <sup>3</sup> 2 *P* <sup>2</sup> Λ*isq*, so it can be controlled through *isq*. The overall control scheme is depicted in Figure 6.

**Figure 6.** Control of the PMSG.

#### 4.2.2. Control for Proportional Load Sharing

As explained in Section 3, when the microgrid is islanded, if the battery is full, a full charge and the DG units' generated power exceeds the total load power demand, and the excessive produced energy starts to flow to the microgrid DC bus capacitors, making the DC voltage rise. A scheme can again be applied to reduce the PMSG output power proportionally to the DC bus voltage increase per unit of the maximum allowed voltage increment.

Based on the wind generation efficiency characteristic shown in Figure 7, it can be observed that if the tip speed ratio *λ* is increased from the optimum value for MPP (*λ opt* = 6.9), the efficiency factor *C<sup>p</sup>* decreases from the maximum value (0.42), becoming zero when *λ* ∼= 15.5. Therefore, an adjustment in the value of the tip speed ratio can be made, according to Figure 8, to force a reduction of generated power in proportion (approximately) to the increase in the DC bus voltage increase.

**Figure 7.** Wind turbine typical efficiency characteristic, considering the fixed blades' pitch angle.

**Figure 8.** Wind generation schemes for MPPT or adaptation to load power demand.

It is essential to mention that, similarly to the PV power reduction method proposed, wind power is reduced from the level corresponding to the MPP condition in a linear relation with the microgrid DC bus voltage increase per unit of the maximum DC bus voltage variation. Thus, the total load power-sharing among all DG units will be approximately proportional to their respective instantaneous MPPs.

#### *4.3. Control of the Microgrid Inverter*

The main objective of the microgrid inverter is to transfer to the microgrid AC bus the excess of power produced by the DG units connected to the microgrid DC bus (or absorb the deficit). Neglecting the system losses, if the net power *PINV*\_*in* injected into the microgrid DC bus (produced by DG units minus the power consumed by the loads and battery being charged) is higher (lower) than the active power *PINV*\_*out* delivered to the microgrid AC bus by the inverter, the difference is absorbed (delivered) by the microgrid DC bus capacitors, and the microgrid DC bus voltage rises (falls). Therefore, the control

objective of the microgrid inverter can be accomplished by regulating its DC bus voltage. The effect of *PINV*\_*in* and *PINV*\_*out* on the microgrid DC bus voltage can be written as:

$$P\_{\rm INV\\_in} - P\_{\rm INV\\_out} = \frac{d}{dt} \left(\frac{1}{2}\mathcal{C}V\_{\rm DC}^2\right) = \frac{1}{2}\mathcal{C}\frac{d}{dt}\left(V\_{\rm DC}^2\right) \tag{15}$$

According to (15), *PINV*\_*out* can be used to control *V* 2 *DC* (and therefore to control *VDC*). Due to the linear relationship between (*PINV*\_*in* − *PINV*\_*out*) and *V* 2 *DC*, and since the reference voltage *V* ∗ *DC* is constant in a steady-state condition, then a proportional–integral (PI) controller is adopted, taking the error *V* 2∗ *DC* − *V* 2 *DC* as input and the negative of the inverter output power −*P* ∗ *INV*\_*out* as the output control action. *PINV*\_*in* is considered a disturbance to be rejected by the controller.

As a second goal, the reactive power delivered by the microgrid AC side should be maintained equal to zero (*Q*∗ *INV*\_*out* = 0) to avoid unnecessary reactive current components flowing through the inverter.

The above shows that it is necessary to regulate the inverter output of active and reactive powers to achieve the objectives. According to the instantaneous power theory [12], the active and reactive powers delivered to the AC bus of the microgrid can be written, in space-vector *αβ* Clarke coordinates, in terms of the inverter AC output current vector → *i* and microgrid AC bus voltage vector <sup>→</sup> *v BUS* as follows:

$$\stackrel{\rightarrow}{s}\_{INV} = p\_{INV} + jq\_{INV} = \stackrel{\rightarrow}{v}\_{BUS} \stackrel{\rightarrow}{i}' = \left(v\_{BUSa} + jv\_{BUS\beta}\right)\left(i\_{\alpha} - ji\_{\beta}\right). \tag{16}$$

From (16), the inverter output current components necessary to impose the reference values of active and reactive powers are

$$
\begin{bmatrix} \dot{i}\_a^\* \\ \dot{i}\_\beta^\* \end{bmatrix} = \frac{1}{\begin{vmatrix} \stackrel{\circ}{\upsilon}\_{BUS} \\ \stackrel{\circ}{\upsilon}\_{BUS} \end{vmatrix}} \begin{bmatrix} \upsilon\_{BUS\alpha} & \upsilon\_{BUS\beta} \\ \upsilon\_{BUS\beta} & -\upsilon\_{BUS\alpha} \end{bmatrix} \begin{bmatrix} P\_{INV\\_out}^\* \\ 0 \end{bmatrix} \tag{17}
$$

Alternatively, using complex space vector notation:

$$
\stackrel{\rightarrow}{i}^\* = P\_{INV\\_out}^\* \frac{\stackrel{\rightarrow}{\vec{v}\_B}\_{BUS}}{\left| \stackrel{\rightarrow}{\vec{v}\_B}\_{BUS} \right|^2}. \tag{18}
$$

Assuming the reference active power *P* ∗ *INV*\_*out* is constant (or slowly varying), if the voltages at the microgrid AC bus are not balanced or contain harmonic components, the currents calculated through (18) would contain the same levels of unbalance and harmonic contamination. For this reason, it is generally preferred to calculate the reference current vector using the positive-sequence fundamental-frequency (FFPS) vector component of → *v BUS*:

$$
\stackrel{\rightarrow}{i}^\* = P\_{INV\\_out}^\* \frac{\stackrel{\rightarrow}{\upsilon}^{+1}\_{BUS}}{\left| \stackrel{\rightarrow}{\upsilon}^{+1}\_{BUS} \right|^2}. \tag{19}
$$

Using (19), the active and reactive powers delivered to the microgrid AC bus might contain oscillations, but their mean values correspond to the reference components. Furthermore, the inverter output phase currents will be balanced and sinusoidal. Finally, to ensure the inverter's desired behavior, inner control loops regulate the inverter *α* and *β* current components. The inverter is connected to the microgrid AC bus through an inductive filter, modeled as an LR circuit where L and R are the filter inductance and

internal resistance. The following vector equation expresses the relationship between the filter voltages and current in the stationary *αβ* reference frame:

$$
\stackrel{\rightarrow}{\vec{v}}\_{INV} - \stackrel{\rightarrow}{\vec{v}}\_{BUS} = R\stackrel{\rightarrow}{\mathbf{i}} + L\frac{d\stackrel{\rightarrow}{i}}{dt}.\tag{20}
$$

Thus, the inverter output voltage vector <sup>→</sup> *v INV* can be used to regulate the current vector → i . In this case, the microgrid AC voltage vector <sup>→</sup> *v BUS* is a known disturbance (since it is necessarily measured for computing the current reference), allowing a feedforward compensation.

In [13,14], several current controllers used in grid-connected converters are described and compared. Advanced control techniques are necessary whenever the voltage at the grid point of common coupling (PCC) contains harmonic components and/or unbalance once the controller is required to reject these disturbances. In these cases, the most popular control strategies recommend using PI controllers in multiple DQ reference frames [15], second-order PR controllers in parallel [16], PI controllers with resonant controllers in a rotating DQ reference frame [17], stationary-frame controllers using the space-vector Fourier transform [18], or repetitive controllers (in a real domain [19] or a complex domain [20–22]). This work used a proportional action and some resonant controllers in parallel due to

the acceptable performance and relatively simple implementation.

The external loop of the inverter control scheme is shown in Figure 9. If the internal controller is designed to be considerably faster than the external one, then the external controller can be designed neglecting the internal dynamics, i.e., considering → *i* / → *i* ∗ = 1.

**Figure 9.** Simplified block diagram of the microgrid inverter's external control loop.

The block diagram of the described internal control loop (*α* and *β* current controllers) is presented in Figure 10.

**Figure 10.** Simplified block diagram of the internal current control loop of the microgrid inverter.

The internal current control loop must have a high bandwidth and high gains at the fundamental frequency and the harmonic frequencies corresponding to typical harmonic disturbances of the grid. These specifications can be met using proportional action and second-order sinusoidal integrators (SSIs) as resonant controllers in parallel. The selection of the appropriate proportional gain enables the imposition of the desired dynamic response to the system, while the resonant terms, implemented in parallel, guarantee a low steadystate error.

A resonant term is included at the fundamental frequency to track the FFPS reference current with a low steady-state error. Resonant terms are also implemented in the low-order harmonic frequencies (3rd, 5th, 7th, and 9th), considering that these are usual components of disturbances in the voltage. In this context, when selecting the maximum harmonic component h = 9, the minimum 0 dB crossover frequency (fPM) required is 540 Hz (an fPM around 1500 Hz was adopted).

In addition to the phase margin (PM) and the gain margin (GM), the sensitivity index *η*, which defines the minimum distance from the Nyquist diagram to the critical point (−1, 0), is also verified. As a design criterion, we considered GM ≥ 3 dB, *η* ≥ 0.3, and PM ≥ 30◦ . The proportional controller gain needed to ensure the desired crossover frequency is *k<sup>p</sup>* = 0.09899. Then, the SSI-based resonant controllers are added in parallel and tuned with the harmonic components h = 6k ± 1 until h = 9. The same resonant gain of *k<sup>i</sup>* = 1 was chosen for all resonant units, and it was observed that the desired phase margin was maintained. Concerning the inevitable delay in digital implementation, a computational delay of two samples was considered. The Bode and Nyquist diagrams of the resulting current control loop are presented in Figure 11a,b, respectively, demonstrating that the design criteria were met.

**Figure 11.** The Bode and Nyquist diagrams of the resulting current control loop. (**a**) Nyquist diagrams of the open-loop system with computational time delay compensation. (**b**) Bode diagrams of the open-loop system with computational time delay compensation. Results were obtained for the proportional controller (blue) and P-SSI in parallel (green).

After defining the internal current control loop gains, the much slower external DC bus voltage loop tuning was considered. The 0 dB crossover frequency of the external loop was chosen so that

$$f\_{cv} < \frac{f\_{ci}}{10} \tag{21}$$

where *fcv* and *fci* represent the 0 dB crossover frequencies of the external and internal control loops, respectively. This design criterion is selected so that the DC bus voltage regulation does not compromise the dynamics of the current control. Therefore, for *fci* = 1500 Hz, *fcv* < 150 Hz must be chosen. For practical reasons, in order to avoid the too frequent saturation of the controller output, *fcv* = 6 Hz was adopted. Furthermore, the indices GM ≥ 3 dB and PM ≥ 30◦ were considered as design specifications. The frequency response of the open-loop transfer function for *k<sup>i</sup>* = 1.0562 and proportional controller gain

*k<sup>p</sup>* = 0.084, which is the value that results in the desired crossover frequency *fcv* = 6 Hz, is shown in Figure 12. The compensated system has a gain margin GM = 65 dB and phase margin PM = 71.5◦ , meeting the design specifications.

**Figure 12.** Bode diagrams of the DC bus voltage open-loop control system.

#### *4.4. Synchronization with the Electrical Network*

In many practical situations, the distribution network exhibits distorted and unbalanced voltages. In these cases, synchronization with the network strongly influences the performance of the entire control scheme of the power converter connected to the network. Determining the correct value of the fundamental-frequency positive-sequence (FFPS) component of the microgrid bus voltage vector is essential for reasonable control of the inverter. This information is obtained through a phase-locked loop (PLL) algorithm.

Several PLL schemes for three-phase systems have been proposed in recent years, the synchronous reference frame PLL (SRF-PLL) being the most popular [23,24]. However, due to its poor performance when the voltages contain harmonics and/or unbalance, preferable alternatives have been presented [25]. In this paper, the PLL scheme that uses the generalized delayed signal cancelation (GDSC) method as a pre-filter to eliminate the effects of unbalances and harmonic components eventually present in the microgrid AC bus voltage was used [26,27]. The choice was based on the fast and accurate determination of the FFPS component of a three-phase signal, even under severe unbalanced and distorted conditions.

#### **5. Operational Scenarios Used for the Validation of the Proposed System**

First, the power flow balance is assumed, i.e., regardless of the state of the distributed generation sources, battery, and load:

$$\mathbf{P}\_{\rm PV} + \mathbf{P}\_{\rm WT} \pm \mathbf{P}\_{\rm ESS} - \mathbf{P}\_{\rm load} \pm \mathbf{P}\_{\rm grid} = \mathbf{0} \tag{22}$$

The power generated photovoltaic panel = PPV; the power generated by the mini-wind generator = PWT, the power injected or absorbed by the storage system = PESS, the total load power = Pload, and the power injected into or absorbed from the electrical network of the university = Pgrid.

Depending on the microgrid operational condition (interconnected or isolated), the instantaneous generated power, the total load power, and battery state-of-charge, the converters must behave according to the modes of operation described.

For interconnected conditions, it is assumed that the battery does not inject any power into the electrical network; in addition, the PV and WT are maintained at their maximum power points.

• Interconnected mode, DG > Pload:

$$\mathbf{P}\_{\rm PV} + \mathbf{P}\_{\rm WT} = \mathbf{P}\_{\rm load} + \mathbf{P}\_{\rm grid}; \mathbf{P}\_{\rm ESS} = \mathbf{0};\tag{23}$$

• Interconnected mode, DG = Pload:

$$\mathbf{P\_{PV}} + \mathbf{P\_{WT}} = \mathbf{P\_{load}} \mathbf{;} \mathbf{P\_{grid}} = \mathbf{0} \mathbf{;} \mathbf{P\_{ESS}} = \mathbf{0} \mathbf{;} \tag{24}$$

• Interconnected mode, DG < Pload:

$$\mathbf{P}\_{\rm PV} + \mathbf{P}\_{\rm WT} + \mathbf{P}\_{\rm grid} = \mathbf{P}\_{\rm load} \text{; } \mathbf{P}\_{\rm ESS} = \mathbf{0}.\tag{25}$$

For isolated conditions, no power can be transferred to or absorbed from the electricity grid. Therefore, if the balance between generation and load demand does not occur, the battery must supply or absorb the difference. However, the DG units must reduce their generated power to meet the load demand if the battery is completely charged. Thus, the following conditions must be verified:

• Isolated mode, DG > Pload, battery not fully charged:

$$\mathbf{P\_{PV}} + \mathbf{P\_{WT}} - \mathbf{P\_{ESS}} = \mathbf{P\_{load}} \text{; } \mathbf{P\_{grid}} = \mathbf{0} \text{, PV, and wind units under MPPT;} \tag{26}$$

• Isolated mode, DG > Pload, battery fully charged:

$$\text{P}\_{\text{PV}} + \text{P}\_{\text{WT}} = \text{P}\_{\text{load}}; \text{P}\_{\text{grid}} = 0; \text{P}\_{\text{ESS}} = 0, \text{PV}, \text{and wind units under reduced power mode}; \quad (27)$$

• Isolated mode, DG = Pload:

$$\mathbf{P}\_{\rm PV} + \mathbf{P}\_{\rm WT} = \mathbf{P}\_{\rm load}; \mathbf{P}\_{\rm grid} = \mathbf{0}; \mathbf{P}\_{\rm ESS} = \mathbf{0}, \text{PV}\_{\rm r} \text{and wind units under MPPT}; \tag{28}$$

• Isolated mode, DG < Pload:

$$\mathbf{P\_{PV} + P\_{WT} + P\_{ESS} = P\_{load}} \mathbf{\hat{P}\_{grid}} = \mathbf{0} \text{, PV, and wind units under MPPT.} \tag{29}$$

#### **6. Results**

A simulation was performed in which load steps were applied for evaluating the operation of the microgrid under all seven of the operational scenarios described.

The behavior of the microgrid DC bus voltage for the different control strategies and operating modes is shown in Figure 13. It can be seen how the DC bus voltage tends to the reference value, regardless of the interconnected or isolated condition of the microgrid, except when the microgrid is isolated, the battery is fully charged, and the load demanded power is lower than the DG available power. In this case, the power generated by each DG unit is reduced in approximate proportion to the increase in the DC bus voltage, as expected. As shown in Figure 14, a relatively smooth transition is achieved for each operating mode of the microgrid.

**Figure 13.** Voltage DC bus microgrid.

**Figure 14.** Relation of electrical power between the control strategies and modes of operation of the microgrid with PPV, PWT, PESS, Pload, and Pgrid.

Initially, the microgrid is connected to the main power system, and three possible situations are considered. (i) DG > Pload: in this case, since the main microgrid DC/AC converter transfers to the grid the excess of generated power, the DC bus voltage is controlled to remain at the rated value, as shown in Figure 13. Additionally, the PV and wind generation units operate under MPPT mode, producing around 480 W and 1900 W, respectively, as observed in Figure 14. The battery is considered to be fully charged and therefore does not absorb any power. (ii) DG = Pload: at t = 3 s, the load power abruptly increases, becoming approximately equal to the total power produced by the DG units. Due to the net power (Pnet = DG-Pload) reduction, the DC bus voltage starts to decrease, but the DC bus voltage controller asks for the DC/AC converter power (Pgrid) to rise until the DC bus voltage returns to the rated value, which occurs when Pgrid = Pnet = 0. (iii) At t = 6 s, the total load power increases again, making DG < Pload. In this case, the DC bus voltage controller is forced to absorb active power from the AC mains (Pgrid < 0) to make VDC = 250 V. The microgrid is connected to the main power system and the DG units generated powers are constant since they remain under MPPT operation. (iv) At t = 9 s, the microgrid interconnection circuit breaker is turned off, making Pgrid = 0. At this moment, the battery interface converter assumes the responsibility for the DC bus voltage control and starts to deliver the necessary net power to regulate the DC bus voltage at its nominal value. (v) At t = 15 s, the total load power decreases so that DG = Pload again. Once the battery is delivering some power to the DC bus voltage, its value tends to increase right after t = 15 s, but the DC bus voltage controller commands a decrease in the power delivered by the battery. (vi) At t = 18 s, an additional decrease in the total load power makes DG > Pload again, and the DC bus voltage experiences a transient increase, but the battery converter again regulates this voltage and absorbs the excess power. However, at some instant after t = 20 s (vii), the battery becomes wholly charged (SoC = max), and the battery cannot absorb any more power. Since DG > Pload, the net power is delivered to the DC bus, the voltage of which naturally starts to increase. During phases (iv), (v), and (vi), the DG units operate under MPPT, but the energy management strategy commands off-MPPT behavior when SoC = max. The power generated by each DG unit decreases in inverse proportion to the DC bus voltage until the total generated power becomes equal to the power demanded by the microgrid loads. As can be seen in Figures 13 and 14, the proposed energy management strategy, together with the proposed DG units' control schemes, operate as theoretically predicted.

#### **7. Conclusions**

The control of a microgrid makes it possible to take full advantage of the distributed generation resources, and this, in turn, allows it to be used to satisfy the energy demand of low voltage electrical installations. These control strategies are necessary to achieve smooth transitions and interaction between the electrical network and the micro-network for each of the operating modes. This paper simulates and evaluates a control strategy that allows power transitions between the microgrid interconnected to the university campus and all possible operating scenarios. The converters have a low number of switches but allow adequate control of each element of the microgrid so that the GD units operate under maximum power point tracking (MPPT) mode or promotes the balance of the generated and load powers.

As observed in the simulations, the correct control strategies for all the elements that make up the distributed generation system, such as for PV MPPT control or load powersharing, PMSG-based wind turbine MPPT control or load power sharing, the control of the microgrid inverter and the synchronization algorithm demonstrated that the proposed strategy is adequate to integrate electrical microgrids in the low voltage network. This type of technology allows a change in the current electrical energy distribution system, making end-users important actors in regulating, controlling, and decongesting electrical networks.

Currently, work is being done on the implementation phase of the proposed control system in the microgrid that is being developed at the Technological Institute of Santo Domingo as part of a project seeking to integrate electrical microgrids in distributed low voltage electrical networks. Future work will be dedicated to implementing the energy management and control system (EMCS) to allow the management of the microgrid in a testing center that will help to analyze the behavior of the control system proposed.

**Author Contributions:** Conceptualization, M.A.-M., L.L.-V. and F.N.; methodology, M.A.-M., L.L.-V. and F.N.; software, M.A.-M. and F.N.; validation, F.S., F.N. and V.A.G.; formal analysis, M.A.-M., L.L.-V. and F.N.; investigation, M.A.-M. and L.L.-V.; writing—original draft preparation, M.A.-M.; writing—review and editing, F.S., F.N., D.M.-H. and V.A.G.; visualization, M.A.-M. and D.M.-H.; supervision, L.L.-V. and F.N.; funding acquisition, M.A.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by FONDOCYT Grant No. 2018-2019-3C1-160 (055-2019 INTEC) in the Dominican Republic.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors gratefully acknowledge the National Background of Innovation and Scientific and Technological Development of the Dominican Republic (FONDOCYT) for the financial support for this research.

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

#### **References**


## *Article* **Linear Programming Coordination for Overcurrent Relay in Electrical Distribution Systems with Distributed Generation**

**Daniel Alcala-Gonzalez 1 , Eva M. García del Toro 1 , M. Isabel Más-López 2 , Sara García-Salgado 1 and Santiago Pindado 3, \***


**Abstract:** Electric power distribution networks are generally radial in nature, with unidirectional power flows transmitted from the highest voltage levels to the consumption levels. The protection system in these distribution networks is relatively simple and consists mainly of fuses, reclosers (RC) and overcurrent relays (OCRs). The installation of distributed generation (DG) in a network causes coordination problems between these devices, because the power flows are no longer unidirectional and can flow upstream to the substation. For this reason, the work proposed here analyzes the most significant impacts that DG has on the protection devices and proposes an adjustment method for the OCRs based on linear programming (LP) techniques with the aim of improving their response time to the different faults that may occur in the main feeder of the network. The distribution system selected for the study is the IEEE 34 bus system using DIgSILENT 14.1 software for its modeling and Matlab for the adjustment of the overcurrent devices. Results indicate that better coordination between protection devices are achieved if LP is used.

**Keywords:** adaptive protection; distributed power generation; power distribution; power system protection

#### **1. Introduction**

The current philosophy in the planning, management and control of a radial power distribution network is based on the assumption of the existence of unidirectional power flows, which is transmitted from the highest transport voltage levels to the distribution levels. We can assume that short-circuit currents behave similarly. These assumptions allow for the implementation of relatively simple and inexpensive protection schemes with which a selective operation of the protection system is achieved [1,2]. According to the principles of selectivity [3], only the protection device closest to the fault should operate to clear the fault, leaving the rest of the network energized.

The installation of DGs at medium and low voltage levels changes this fundamental basis. Power flows and short-circuit currents can now have different upstream directions and values [4]. As a consequence, the initial schemes implemented (e.g., main feeder protection) may no longer work or be less effective [1,5].

With the presence of DG the distribution system becomes more active, i.e., both load and generation significantly affect the state of the network [6]. The effects of DG on distribution network devices have been discussed in different literature [4,7]. To reduce them, several "mitigation methods" have been proposed, thus, several authors recommend acting directly on the DG [3,8], e.g., disconnecting it just before fault detection or limiting its

**Citation:** Alcala-Gonzalez, D.; García del Toro, E.M.; Más-López, M.I.; García-Salgado, S.; Pindado, S. Linear Programming Coordination for Overcurrent Relay in Electrical Distribution Systems with Distributed Generation. *Appl. Sci.* **2022**, *12*, 4279. https://doi.org/ 10.3390/app12094279

Academic Editors: Luis Hernández-Callejo, Maria del Carmen Alonso García and Sara Gallardo Saavedra

Received: 27 March 2022 Accepted: 19 April 2022 Published: 23 April 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/).

power and therefore its effect on the protection system [9], others recommend reconfiguring the network topology [10,11], or the use of fault current limiters (FLCs) [12,13].

In recent years, mitigation methods based on coordination algorithms have been developed using linear programming (LP), non-linear programming (NPL), or genetic algorithms (GA) [14–17] that either modify the coordination interval time (CTI) between devices, determine the optimal location of the DG [15], or in combination, without compromising the protection system.

Although these methods can be effective, they have some disadvantages and limitations. Therefore, the disconnection of the DG immediately before fault detection proposed in [3,18], can cause asynchronous reconnections and cause severe damage to both the DG and the distribution network, especially at high penetration levels [3,10]. The limitation of the DG capacity proposed in [9] is also undesirable, since it also limits the penetration level and therefore the advantages of its installation close to the consumption points, (loss reduction, improvement of supply quality, etc.). The network reconfiguration suggested in [10,11] is costly and sometimes unfeasible. Finally, the use of current limiters [12,14] entails an additional installation cost that can be prohibitive.

Regarding the use of NLP methods for the determination of the coordination index mentioned in [15] can be very complex, particularly if the number of protection devices and DG's is high, since the calculation time required to determine their operating parameters will also be high [19–21], on the other hand, not always the optimal location point of DG is viability from the perspective of network operation, e.g., availability of power evacuation at the point of common connection (PCC).

The investigations described in [14,16,17], use GA-based optimization algorithms to determine optimal coordination between OCR and distance relays [14]. In [16] a hybrid GA and NLP algorithm is developed to choose the optimal value of DIAL or time multiplier setting (TMS) of all OCR, the use of a continuous genetic algorithm (GAC) faster than GA for the same purpose is developed in [17]. However, all the optimization algorithms described do not consider DG.

In the present work:


As far as the authors know, it seems that there is a lack of similar studies in the available literature.

#### **2. Effect Caused by the DG on Protective Devices**

#### *2.1. Loss of Sensitivity*

The problem related to the loss of sensitivity to main feeder protection is often referred to as protection blinding [8]. Installation of DGs in the distribution system can reduce the value of the short-circuit current detected by the protection of the main substation and affect the response time of the circuit breaker, which will depend, to a large extent, on the size and location of the DG within the distribution network.

The loss of sensitivity of the protection device can be analyzed as a function of where the short circuit is located in relation to the location of the DG, distinguishing between short circuits located upstream of the DG and short circuits located downstream of it. Figure 1 shows the topology of a radial distribution network where S.E is the substation, R is the main feeder protection, RC is the automatic recloser, and F is a protection fuse for the branch line.

**Figure 1.** Sketch of a short circuit upstream of DG.

For faults located upstream (between the substation—S.E, and the DG), the contribution current of the S.E to the short-circuit is independent of the size of the DG. A three-phase short-circuit is also represented upstream of the DG, for different penetration levels: 17%, 33%, and 50%. In the operating characteristic curves of the time overcurrent (t/i) devices shown in Figure 2, this effect is observed, as the current detected by the relay remains constant at 1.59 kA, regardless of the of penetration level of the DG.

**Figure 2.** Short circuit upstream of DG. See also Figure 1.

In contrast, for faults located downstream of the substation and the DG (as shown in Figure 3), it is observed that the current contribution of the DG to the short circuit causes the current measured by the main relay, R, of the S.E to decrease from 1.21 with 17% penetration to 1.14 kA with 50% penetration level of the DG, which represents a 6.1% loss in the sensitivity of the protection device. This effect will cause a delay in the operating time of the protection device, as shown in Figure 4. The extreme case of this effect occurs when the main relay does not detect the fault current, as shown in Figure 5. The maximum value at which the relay stops detecting the fault current appears when the DG is installed close to the substation and under fault conditions at the main feeder [12].

**Figure 3.** Sketch of a short circuit downstream of DG.

**Figure 4.** Short circuit downstream of DG. See also Figure 3.

**Figure 5.** Short loss of sensitivity substation relay. Loss of sensitivity.

#### *2.2. Loss of Coordination*

Under normal operating conditions, protection devices are coordinated in such a way that, in the event of a fault in the network, the main protection acts before backup protection. A large majority of faults occurring in the network are of a temporary nature, the purpose of the recloser, RC, is to try to clear these faults and on the one hand avoid an unnecessary interruption of the electrical service, and on the other hand to safeguard the protection fuse of the branch line under its supervision. Depending on the initial settings of the coordination schemes, and taking into account the size, location, and type of DG installed, it may happen that a temporary fault in a branch line causes the loss of coordination between fuse and recloser, affecting their coordination time due to the additional current with which the DG contributes to the short-circuit. To observe this impact, the situation of temporary fault, ICC, in a branch line (lateral) is shown in Figure 6.

**Figure 6.** Short circuit in branch line.

The characteristic t/i trip curves are represented in Figure 7 where the impact called (fuse nuisance blowing) is shown. We can observe how a temporary fault in a branch line causes a loss of coordination between the fuse and the recloser. Thus, the fuse detects a fault current of 0.25 kA and trips in 0.15 s, while the recloser is crossed by a current of 0.650 kA, which would cause it to trip in 0.32 s. This time is higher than the fuse tripping.

**Figure 7.** Failure coordination between recloser, RC, and fuse, F, see also Figure 6.

#### **3. Coordination with Linear Programming**

OCRs are currently the most widely used protection devices [13,22], this device generally being used as backup protection. However, in some situations, it may be the only protection. The OCR should act as primary protection by removing faults that occur in the area under its supervision. Only in the event that the primary protection fails to operate should the backup protection initiate tripping. A typical electrical distribution network may consist of hundreds of equipment protection relays. Each relay in the system must be coordinated with the adjacent equipment protection relay. If backup protection is not well coordinated, coordination failure may occur, and therefore, coordination of backup protections is one of the main concerns in the network protection system [13,23]. Generally, protection coordination can be performed by topology [24,25], by optimization methods [26,27], or by expert methods [28]. Topological analysis is used for relay tuning in multiterminal networks, graph theory, and functional approximation techniques are employed to provide the best solution that does not necessarily have to be the optimal one. In optimization methods, some researchers [26,29] use non-linear programming techniques to determine the optimal relay settings, subject to constraints due to coordination and limits of the relay settings themselves. In reference [30], the big M method is proposed to find the optimal value of the OCR time multiplier setting (TMS) in which the stated values of the plug setting are assumed to be known and fixed.

The problem of OCR coordination in the distribution system with DG presence can thus be defined as an optimization problem with constraints. The objective is to minimize the operating time of the relay closest to the place where the fault has occurred. The constraints imposed are due to limits on the operating time of the relay, coordination criteria, and relay characteristics. In this work, the OCR coordination problem in radial distribution systems is formulated as a linear programming (LP) problem with constraints, where the sum of the operating times of the relays in the system for different minimum fault points:

$$\min \frac{\alpha}{\left(k\_{i,j}\right)^{\beta} - 1} \sum\_{i=1}^{n} \left(\sum\_{j=1}^{m} \text{TMS}\_{i,j} - \text{TMS}\_{i+1,j+1}\right) . \tag{1}$$

where *α* and *β* are OCR shape constants according to Table 1, *ki,j* is the ratio between the short-circuit current, *ICC*, and the relay setting current, *Iar*, regarding relay *i* with fault in section *j*, and TMS is the time multiplier setting.


**Table 1.** Form constants for the exponential equation IEC [31].

The proposed optimization problem is subjected to the following group of constraints:

• Coordination criterion—The protection coordination criterion establishes the minimum time that must elapse between the operation of the primary protection and the operation of the backup protection. The fault is detected simultaneously by the primary protection, PP, and the secondary protection, PS. To avoid erroneous operation, the PS will only have to operate in case the PP fails. If we define *R<sup>i</sup>* as the primary fault protection at a certain point, *j*, and *Ri*+1 as the secondary or backup protection for the same fault. The constraint condition for coordination criteria is:

$$t\_{i+1,j} - t\_{i,j} > \Delta t\_{\prime} \tag{2}$$

where, *ti*+1,*<sup>j</sup>* is the operation time of the back-up, PS, operation for fault at point *j*, *ti,j* is the operation time of PP operation for the same fault, and ∆*t* is the coordination time interval (CTI).

• Relay operating time limits—The time taken by a relay to detect and isolate a fault produced in its zone of influence must be bounded. This is the constraint imposed by the operating time of the relays:

$$t\_{i,min} \le t\_{i,j} \le t\_{i,max} \tag{3}$$

where, *ti*,min and *ti*,max are the minimum and maximum operating times of relay *i* (at any point).

• Time multiplier setting, TMS, limits—The operating time of a relay is directly proportional to the TMS. Therefore:

$$\text{TMS}\_{i,\text{min}} \le \text{TMS}\_{i,\text{j}} \le \text{TMS}\_{i,\text{max}} \tag{4}$$

where TMS*i*,min is the minimum TMS value for relay *i* and, TMS*i*,max is maximum value of TMS for relay *i*. The TMS values usually taken are 0.025–1.2, respectively [32].

• Relay operating characteristics—To extend the work region of the overcurrent relay specified by the standard [33], the relay's parameters are optimized considering the maximum value of the relay and the minimum value of the relay, as expressed respectively. We consider all relays with the same characteristics:

$$t\_{op} = \frac{\alpha}{\text{PSM}^{\beta} - 1} \text{TMS}\_{\text{\textdegree}} \tag{5}$$

where *top* is the relay operating time and PSM is the plug setting multiplier (*ki,j* in Equation (1)). For inverse characteristic curves it is usually *α* = 0.02 and *β* = 0.14 [32,33].

The relay setting current, *Iar*, is determined from the system requirements. Thus, the above equation can be rewritten as:

$$t\_{op} = a \times \text{TMS}\_{\prime} \tag{6}$$

where:

$$a = \frac{\mathfrak{a}}{\text{PMS}^{\beta} - 1}.\tag{7}$$

#### **4. Case Study Results and Discussion**

The distribution system selected for the simulations is the one taken from the test feeders of the Distribution System Analysis Subcommittee of the Institute of Electrical and Electronics Engineers (IEEE) [34], shown in Figure 8. The main characteristics of the distribution network on which the short circuits were simulated are the following:


**Figure 8.** Modified IEEE 34-node test feeder system [34].

Three-phase and single-phase short circuits were simulated on the distribution network, at busbars 2, 7, and 21, as shown in Figure 8, and with different DG penetration levels ranging from 33% to 50%. A first distributed generation unit, DG1, was installed at busbar 20, whereas a second one, DG2, was installed at busbar 22.

The DG protection device should be properly coordinated with the other protection devices in the distribution network. When a fault occurs on the line to which the DG unit is connected, the DG unit must be disconnected before the protection device on that line is disconnected. This ensures that the DG does not cause disturbances in the operation of the line protections. Furthermore, when a fault occurs on a branch line, the DG must remain in service and will not disconnect before the fault line does [8]. According to this, the coordination criterion shown in Table 2 is proposed.


**Table 2.** Coordination criteria according fault location.

The transformation ratio of the current transformers used in the test is 400/5, chosen according to the full-load currents of the system and the short-circuit currents detected by each of the protection devices according to the location of the fault (shown in Table 3).


**Table 3.** Short circuit (expressed in amperes, A) detected by the protection devices, for different grounding resistor values, *R<sup>f</sup>* , expressed in ohm.

#### *4.1. Objective Function According to Fault Location*

According to the location of the fault, and based on the protection criteria included in Table 2, three objective functions to minimize can be obtained, each of them will result in a different adjustment, an adaptive adjustment.

The scenario analyzed is the three-phase short-circuit scenario with a DG1 penetration level of 17%, the results obtained for the other scenarios can be found in Table A3 in the Appendix A.

Fault 2: OCR<sup>2</sup> with OCR<sup>1</sup>

$$\text{minz} = \sum\_{i=1}^{m} a\_{i,k} (\text{TMS})\_i = 2.27 \text{ TMS}\_1 + 3.33 \text{ TMS}\_2 \tag{8}$$

Fault 7: Recloser with OCR<sup>2</sup>

$$\text{minz} = \sum\_{i=1}^{m} a\_{i,k} (\text{TMS})\_i = 3.49 \text{ TMS}\_2 + 2.21 \text{ TMS}\_{\text{RC2}} \tag{9}$$

Fault 21: OCR<sup>3</sup> with recloser


$$\begin{array}{l} \text{OCR1 } t\_{\text{op1}} = a\_1 \text{(TMS}\_1) = 2.27 \cdot 0.17 = 0.39 \text{ s} \\ \text{OCR2 } t\_{\text{op2}} = a\_2 \text{(TMS}\_2) = 3.33 \cdot 0.030 = 0.10 \text{ s} \end{array} \tag{18}$$

Fault 7: Recloser with OCR2:

$$\begin{aligned} \text{OCR2 } t\_{\text{op}2} &= a\_2 (\text{TMS}\_2) = 3.49 \cdot 0.038 = 0.13 \text{ s} \\ \text{RECLOSER } t\_{\text{op}\text{RC}} &= a\_{\text{RC}} (\text{TMS}\_{\text{RC}}) = 2.21 \cdot 0.11 = 0.24 \text{ s} \end{aligned} \tag{19}$$

Fault 21: OCR<sup>3</sup> with Recloser

$$\begin{aligned} \text{OCR3 } t\_{\text{op3}} &= a\_3 (\text{TMS}\_3) = 4.57 \cdot 0.065 = 0.30 \text{ s} \\ \text{RECLOSER } t\_{\text{opRC}} &= a\_{\text{RC}} (\text{TMS}\_{\text{RC}}) = 3.23 \cdot 0.025 = 0.1 \text{ s} \end{aligned} \tag{20}$$

The values of constants *a* and TMS and relay operation times *top* for the rest of the scenarios are shown in Appendix A Tables A1 and A2. Table A3, shows a comparison of the operation times between the coordination performed in a classical way and with the calculation of adjustments carried out by using LP. Implementing the results obtained in the simulation software, we obtain the following protection schemes shown below.

The classical coordination plots versus linear programming coordination in the case of three-phase short-circuits in the locations described with 17% DG1, are shown in Figures 9–11. It can be seen that the response of the OCR protection devices is optimized using the method proposed in this work.

**Figure 9.** Short-circuit in F2 with penetration level of DG1 17 %. Classical coordination versus linear programming coordination.

**Figure 10.** Short-circuit in F7 with penetration level of DG1 17%. Classical coordination versus al linear programming coordination.

**Figure 11.** Short-circuit in F21 with penetration level of DG1 17%. Classical coordination versus linear programming coordination.

#### **5. Conclusions**

In this work, a unified protection system that offers double functionality is proposed. The system was proven to be capable of both optimizing the settings of the relays for each network situation.

The performance of the relay has been tested for different DG penetration levels, different types of short-circuit faults (single-phase to ground and three-phase), and different fault zones (faults in the main feeder and the laterals).

The proposed scheme can be implemented in distribution systems using the communication standard IEC61850 [31].

In all the cases studied, the system has proven capable of adapting, in real time, the performance parameters of the relays that protect the faulted line section depending on the operating conditions of the system at the time of the fault.

The proposed protection system has been applied to a distribution system in the presence of DG. The results demonstrate that the protection coordination proposed in this paper is capable of reducing relay operating times by more than 80% compared to classical coordination, which produces a more reliable operation of a network with DG.

**Author Contributions:** Conceptualization, E.M.G.d.T. and M.I.M.-L.; methodology, D.A.-G., E.M.G.d.T., M.I.M.-L., S.G.-S. and S.P.; software, D.A.-G.; validation, S.G.-S., E.M.G.d.T. and M.I.M.-L.; formal analysis, E.M.G.d.T., M.I.M.-L. and D.A.-G.; investigation, D.A.-G., E.M.G.d.T., M.I.M.-L., S.G.-S. and S.P.; data curation, E.M.G.d.T. and M.I.M.-L.; writing—original draft preparation, D.A.-G., E.M.G.d.T., S.G.-S., M.I.M.-L. and S.P.; writing—review and editing, E.M.G.d.T. and M.I.M.-L.; supervision, S.G.-S., M.I.M.-L. and S.P. 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.

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

#### **Appendix A**


**Table A1.** Values of constant relay (*ai,j*) for the different scenarios analyzed.

**Table A2.** Relay TMS values for the different scenarios analyzed.



**Table A2.** *Cont.*

**Table A3.** Operation time PL coordination versus operation time from classical coordination.


C—classical coordination; **PL**—lineal programming coordination.

#### **References**


## *Article* **Effect of Distributed Photovoltaic Generation on Short-Circuit Currents and Fault Detection in Distribution Networks: A Practical Case Study**

**Daniel Alcala-Gonzalez 1 , Eva Maria García del Toro 1 , María Isabel Más-López <sup>2</sup> and Santiago Pindado 3, \***


**Abstract:** The increase in the installation of renewable energy sources in electrical systems has changed the power distribution networks, and a new scenario regarding protection devices has arisen. Distributed generation (DG) might produce artificial delays regarding the performance of protection devices when acting as a result of short-circuits. In this study, the preliminary research results carried out to analyze the effect of renewable energy sources (photovoltaic, wind generation, etc.) on the protection devices of a power grid are described. In order to study this problem in a well-defined scenario, a quite simple distribution network (similar to the ones present in rural areas) was selected. The distribution network was divided into three protection zones so that each of them had DG. In the Institute of Electrical and Electronic Engineers (IEEE) system 13 bus test feeder, the short-circuits with different levels of penetration were performed from 1 MVA to 3 MVA (that represent 25%, 50%, and 75% of the total load in the network). In the simulations carried out, it was observed that the installation of DG in this distribution network produced significant changes in the short-circuit currents, and the inadequate performance of the protection devices and the delay in their operating times (with differences of up to 180% in relation to the case without DG). The latter, that is, the impacts of photovoltaic DG on the reactions of protection devices in a radial distribution network, is the most relevant outcome of this work. These are the first results obtained from a research collaboration framework established by staff from ETSI Civil and the IDR/UPM Institute, to analyze the effect of renewable energy sources (as DG) on the protection devices of a radial distribution network.

**Keywords:** coordination protection; distributed generation; photovoltaic resources; DigSILENT

## **1. Introduction**

The distributed generation (DG) based on renewable energy sources in some distribution networks has caused a major change in the traditional model of power supply. According to the new paradigm of DG, generation of electricity power should come closer to the user, in contradiction to the traditional generation in centralized plants [1,2]. Therefore, it is necessary to implement electrical systems with an infrastructure that allows the energy to be distributed and made available to users in optimal conditions for their use [3].

Although the presence of DG in power networks can be beneficial [4,5], it produces a structural transformation, failing to behave as a classical radial distribution network. Reversed power flows can occur both in steady state and in the presence of faults [6–10]. The impacts that DG may have on these distribution networks will depend, among other aspects, on the type of generation, technology used, installed power, and location in the

**Citation:** Alcala-Gonzalez, D.; García del Toro, E.M.; Más-López, M.I.; Pindado, S. Effect of Distributed Photovoltaic Generation on Short-Circuit Currents and Fault Detection in Distribution Networks: A Practical Case Study. *Appl. Sci.* **2021**, *11*, 405. https://doi.org/ 10.3390/app11010405

Received: 16 November 2020 Accepted: 10 December 2020 Published: 4 January 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/).

network [11–15]. Furthermore, these new fault modes resulting from the introduction of DG, based on renewable energy such as photovoltaic (PV) generation, can affect the reliability of these PV systems and decrease the revenues foreseen [16,17].

Among the different problems arising with the increasing of DG within power networks, the protection relays coordination is one of the most relevant (other problems being harmonic distortion, frequency drop, and stability and reliability of the network [18–22]). Protection in traditional networks was designed for unidirectional power distribution. However, this situation has changed with the greater importance of renewable energy [23], and other ways to reduce power production from the traditional energy sources that have a negative impact on the climate (e.g., the use of electric vehicles as power sources [24–29]). Within the networks with DG, such as the ones with photovoltaic generation located in different points of the grid, the coordination between protection devices might fail [30–38]. Additionally, each kind of renewable energy source presents new fault modes in a network. In [39], there is thorough review of the fault modes that a photovoltaic installation can bring to a distribution network.

The management and control of a radial electric distribution network is based on the assumption of the existence of unidirectional flows of power, which is transmitted from the highest levels of transport voltage to the distribution levels. We can assume that the shortcircuit currents behave in a similar way. These assumptions allow the implementation of relatively simple and economical protection schemes, in order to achieve selective operation of the protection system (according to the principles of selectivity, only the protection device closest to the defect must function to clear the fault, leaving the rest of the network energized). There are numerous literature reviews on the matter of the protection of power grids with DG; the recent works [40–44] being worthy of mention.

The installation of DG in medium and low voltage levels changes the fundamental basis of the aforementioned unidirectional power flow [45]. Both the power flows and the short-circuit currents can now have upstream addresses, these also being of different values than the ones initially foreseen [6]. As mentioned above, the initial schemes applied (that is, the main feeder protection) might start being less effective or even stop working [46]. This problem being caused as the value of the short-circuit current detected by the main protection of the substation may be altered, and therefore affect the response time of the protection devices that will depend, to a great extent, on the size and location of the DG within the distribution network [15,47,48].

As stated, one of the main problems of the presence of DG in distribution networks is the loss of coordination and untimely triggering of protection devices [49], as short-circuit currents can increase as a result of the contribution from those DG. The fault currents suffer variations that can run the fuse in both directions, which means that it can be traversed by currents generated in points downstream and upstream of its location. Therefore, it is necessary to reconfigure the protection devices of these distribution networks. The loss of sensitivity of the protection devices can then be analyzed according to the location of short-circuits in relation to the DG, differentiating among faults located upstream and downstream of the aforementioned DG.

The aim of this work is to describe the preliminary results obtained from a research collaboration framework established by researchers from ETSI Civil and the IDR/UPM Institute to analyze the effect of renewable energy sources (as DG) on the protection devices of a power grid. The aforementioned framework is based on the work already carried out on protection in power grids with distributed generation [50], and the work carried out in analytical modeling on solar panels and batteries [51–60]. In the present work, photovoltaic energy was selected as the example of DG. It should be underlined that this particular source of energy is associated with some specific problems when considered as DG [61–65], and should be analyzed while taking into account its particular performance in relation to temperature and irradiance [66]. Finally, it should be underlined that although there are many works in the available literature on the effect of DG on

distribution grids [6,11,12,14,18,20,24,29,36,67,68], it appears that the effect of DG on the protection coordination of relays curves.

In order to obtain quick results that could lead to general but solid conclusions, a simple distribution network has been analyzed in the present work. This specific power distribution grid was selected as it represents a well-defined scenario in which relevant conclusions can be derived.

As a first work on protection for grids with DG based on renewable energy sources, we chose to analyze a simple distribution network similar to the ones typical of rural areas. With the objective of ensuring the electricity supply in these distribution networks when DG is connected, and maintaining the functionality of the protections, the following aspects were studied:


Finally, it is also fair to say that the present work analyzed the effect of photovoltaic distributed generation in a network, leaving aside other possible renewable sources, such as wind power generation. This was a drawback that will be overcome in future works.

The following paper is organized as follows: in Section 2, the distribution network and different cases studied are described, the results being included in Section 3. Finally, conclusions are summarized in Section 4.

#### **2. Materials and Methodology**

The protection of the distribution networks studied in this work was guaranteed by main protection relays located at the beginning of the line (that is, at the electric substation (S)), and cut-out fuses (F1, F2) at the laterals. The basic topology of this kind of network is shown in Figure 1, where a distributed generation (DG) source has been added. The main feeder is protected by a relay (R), whereas the laterals supplying the loads are protected by means of fuses.

**Figure 1.** Basic topology of the simple power distribution network studied. A distributed generation (DG) source has been added.

Relays are normally equipped with inverse time overcurrent devices, whose performance is defined by the following equation [67]:

$$t(I) = \frac{a}{\left(\frac{I}{I\_{pick-up}}\right)^b - 1} \text{TMS} \tag{1}$$

where *t* is the time of operation of the overcurrent device, *I* is the fault current detected by the device, *Ipick-up* is the adjustment current (relay pick-up curve) of the device, *a* and *b* are the control parameters of the actuation curve, and TMS is the Time Multiplier Setting (expressed in seconds) [69]. The characteristic values of the relays used in the present work are included in Table 1.

**Table 1.** Parameters for the different relay characteristics used in the Institute of Electrical and Electronic Engineers (IEEE)-13 bus test feeder system used in the present work (see Equation (1)). These values were selected according to the IEC 60255-51 standard.


Each relay in a protection system has a set of relay settings that determine the primary and backup protection that the relay will provide. The relay settings for each relay are calculated so that the relay fulfills the primary and backup protection requirements of the network it is protecting. Calculations are based on the maximum load current, the maximum and minimum fault currents, and/or the impedance of feeders that the relay is protecting.

On the other hand, fuses have an inverse current-time characteristic that is usually plotted as a log-log curve, which is better approximated by a second-order polynomial function. However, it should be underlined that the most interesting part of this curve, in terms of its practical applications, can be approximated by a linear expression [70]. Therefore, it can be assumed that the general equation for the characteristic curve of an expulsion fuse can be expressed as [68,71–74]:

$$
\log(t) = m \log(I) + n \tag{2}
$$

where *t* is the actuation time of the ejection fuse, *I* is the current through it, and *m* and *n* are fuse constants to be determined as in [75].

The present work was carried out by using computer simulation of a rural/small grid with 3 possible photovoltaic DG sources, selected as a case study test network. To analyze the effect of the photovoltaic DG on the protection devices of such a distribution network, the IEEE-13 bus test feeder system was used (see Figure 2).

The simulation was performed with DIgSILENT Power Factory software, already used to study the protection of distribution systems with DG [72,76–79]. This distribution network design has been successfully used to study static short-circuit currents [80], the maximum possible photovoltaic power penetration into a network in relation to the demand response [81], voltage regulation strategies in DG [78,82,83], or fault ride-through in power networks related to renewable energy (wind and photovoltaic) [79]. Besides, it should be also said that DIgSILENT Power Factory software is a powerful simulation tool that integrates the following standard methodologies for short-circuit calculations: IEC 60909, IEEE 141/ANSI C37, VDE 0102/0103, G74, and IEC 61363.

The analyzed distribution network was divided into three Zones of Protection (Zones 1, 2, and 3), assigning a specific area for each DG (starting each one of these zones from the beginning of the output feeder from the electrical substation, see Figure 2). The total load of the distribution network was 4 MW, which was considered constant in all simulations. The connection points of the DG were chosen based on their distance from the main substation.

The three different case studies analyzed in the present work are indicated Figure 2 (each of them involves only one photovoltaic DG):

• N1: installation of a photovoltaic DG at bus 646, located at the middle of the distance to the substation (*L* = 2400 m), within Protection Zone 1;


These case studies were analyzed with different DG penetration levels (1, 2, and 3 MW, representing 25%, 50%, and 75% of the power load, respectively). These levels represent normal variations in the behavior of the solar panels in relation to irradiance on the solar cells and their temperature (the use of Maximum Power Point Tracking is supposed). In recent works, we have reflected these variations (both measured and calculated by using implicit and explicit models) [59,60,84].

To the authors' knowledge and after a thorough review of the available literature, the present approach to analyze the performance of a protection scheme based on both relays and cut-out fuses in a well-known power distribution grid, and in relation to the power supplied from photovoltaic DG, represents a novelty.

**Figure 2.** Example of the rural network topology studied in the present work. The position of the three DG included in the topology, N1, N2, and N3 (case studies), is indicated. The different protection zones are also indicated in the figure, together with the simulated short-circuits.

#### **3. Results and Discussion**

The results, as presented for all scenarios, show various ranges of DG capacity which provide different values of fault current. Adding distributed generation to the 20 kV distribution network creates a situation where networks that were designed primarily as tail, radial, or open radial networks become looped networks.

As a result, distributed generation can cause relays in a protection system to underreach or over-reach. This has been illustrated in this paper, by sample calculations and an actual protection review (using DigSILENT PowerFactory software).

In Figure 3, the electric circuit corresponding to Protection Zone 1 with a 1 MW photovoltaic DG source (connected to bus 646) is shown (see also Figure 2). Protection elements in this area include an overcurrent relay at the beginning of the line (bus 632), and a fuse (F645) to protect the line connecting the DG. Only a 1 MW power supply through

bus 646 was analyzed, as it represented the maximum power that was able to be injected by the DG in continuous operation (larger installed powers always tripped fuse F645).

**Figure 3.** Protection Zone 1 electric circuit (see also Figure 2), in which a short-circuit in bus 645 is indicated.

A three-phase short-circuit at bus 645 was simulated; the characteristic curves being plotted in Figure 4. In the graph, it shows that in the case of a fault between the substation and the DG, the presence of the latter represented a decrease in the value of the short-circuit current detected by the relay (R632 in Figure 3), from 536 A to 513 A. Thus, a loss in relay sensitivity (binding), as relay R632 at the main feeder detected a lower value of the fault current, with the triggering time also being delayed from 0.3 s to 0.35 s. These results are also summarized in Figure 5, and Table 2 at the end of this paper.

**Figure 4.** Protection Zone 1 current-time (*I-t*) performance in case of a three-phase short-circuit with 25% DG power (1 MW) at bus 645 (see Figure 3).

**Figure 5.** Triggering time, *t*, of protection relay R632 in Protection Zone 1 (see Figures 2 and 3), in relation to the power supply from DG. The currents detected by the relay when triggering were added to the graph in each case.

**Table 2.** Relay and fuse currents, *I*, and triggering times, *t*, in the protection devices from the defined Protection Zones 1, 2, and 3 (see Figure 2, Figure 3, Figure 6, and Figure 10), in relation to the power supplied by the DG (1, 2, and 3 MW, representing 25%, 50%, and 75% of the photovoltaic installed power).


The Protection Zone 2 analyzed scenario (Figure 2) involved a 1, 2, and 3 MW DG power supply. A protection relay R632 was installed at the substation (see Figure 6), and two faults were studied, at points F1 (downstream from the DG source) and F2 (upstream from the DG source). In case of a short-circuit at point F1, the current detected by the relay decreased from 488 A (without DG) to 455 A (with DG), with a triggering delay of 13 ms (see in Figure 7 the case corresponding to 50% DG power at bus 633). Similarly, in case of a short-circuit at point F2, the currents decreased from 498 A to 467 A, with a triggering delay of 10 ms (see in Figure 8 the case corresponding to 50% DG power at bus 633). The results are summarized in Figure 9 and Table 2.

As performed in relation to the fault analysis carried out in Protection Zone 2, three different DG power supply scenarios (1, 2, and 3 MW) were studied in relation to Protection Zone 3 (Figure 2). The protection devices installed were two overcurrent relays (R671 and R632) within the main distribution line, and a fuse 684 DG protection, see Figure 10. Two three-phase short-circuits were simulated within the corresponding Protection Zone 3: at point F1, located downstream from the DG, and at point F2, upstream from the DG. The results are respectively included in the graphs in Figures 11 and 12. In case of short-circuits

at F1, the DG installation implies a decrease of the current detected by the relay, from 1113 A to 420 A, which represents a 62% decrease of the relay sensitivity and a delay of 550 ms in relation to the triggering time (see Figure 11). Similarly, in the case of short-circuit at point F2 (see Figure 12), the connection of the DG source introduced a decrease of the detected current from 1400 A to 490 A, and a delay of the triggering time, from 300 ms to 800 ms. The results are summarized in Figure 13 and Table 2.

**Figure 6.** Protection Zone 2 electric circuit (see also Figure 2), in which short-circuits in points F1 and F2 are indicated.

**Figure 7.** Protection Zone 2 current-time (*I-t*) performance of in case of a three-phase short-circuit at point F1, with 50% DG power (2 MW) at bus 633 (see Figure 6).

**Figure 8.** Protection Zone 2 current-time (*I-t*) performance in case of a three-phase short-circuit at point F2, with 50% DG power (2 MW) at bus 633 (see Figure 6).

**Figure 9.** Triggering time, *t*, of the main feeder protection relay R632 in Protection Zone 2 (see Figures 2 and 6) for fault F1 downstream from the DG and fault F2 upstream from the DG, in relation to DG power supply. The currents detected by the relay when triggering were added to the graph in each case.

**Figure 10.** Protection Zone 3 electric circuit (see also Figure 2), in which short-circuits in points F1 and F2 are indicated.

**Figure 11.** Protection Zone 3 current-time (*I-t*) performance in case of a three-phase short-circuit at point F1, with 25% DG power (1 MW) at bus 652 (see Figure 10).

**Figure 12.** Protection Zone 3 current-time (*I-t*) performance in case of a three-phase short-circuit at point F2, with 25% DG power (1 MW) at bus 652 (see Figure 10).

**Figure 13.** Triggering time, *t*, of the main feeder protection relay R632 in Protection Zone 3 (see Figures 2 and 10) for fault F1 upstream the DG and fault F2 downstream the DG, in relation to DG power supply. The currents detected by the relay when triggering were added to the graph in each case.

In Table 2, the results obtained in relation to each of the analyzed scenarios are included. The fault currents and the protection relay tripping/triggering times are shown with regard to:


#### **4. Conclusions**

A coordination study was presented and applied to the IEEE 13-node test feeder to evaluate the effect of the PV DG penetration on the protection devices coordination. Different cases were studied by changing DG penetration levels and locations for each possible fault location. A protection coordination assessment was carried out by analyzing the location of the faults (devices see fault currents for downstream and upstream faults).

The most relevant conclusions from this work are:


Additionally, the present protection procedures need to be revised, as the presence of DG sources might represent a quite serious effect on the current and triggering time reaction of the protection relays. According to the results and depending on the distance to the substation, the tripping time with DG was increased up to 167% (for 25% PV DG penetration) and 180% (for 50% and 75% PV DG penetration).

The results from this work can be reasonably extrapolated to grid networks, bearing in mind the current direction effects on the protection relays, and the DG penetration.

Future works of the research team that carried out the present work are the development of simplified techniques to locate faults in distribution networks with DG, new methods to select the appropriate fuses in these networks (to avoid unexpected fuse tripping), and a procedure to successfully reprogram overcurrent relays.

Finally, it should also be noted that the present work will be followed by others that will include research on the dynamic modeling of power supplies from DG in distribution networks, and simplified methodologies on fault locations in grids with DG.

**Author Contributions:** Conceptualization, D.A.-G. and S.P.; methodology, D.A.-G., E.M.G.d.T., M.I.M.-L. and S.P.; software, D.A.-G.; validation, D.A.-G., E.M.G.d.T. and M.I.M.-L.; formal analysis, D.A.-G.; investigation, D.A.-G., E.M.G.d.T., M.I.M.-L. and S.P.; data curation, D.A.-G.; writing original draft preparation, D.A.-G. and S.P.; writing—review and editing, D.A.-G. and S.P.; supervision, D.A.-G. and S.P. All authors have read and agreed to the published version of the manuscript.

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

**Acknowledgments:** The authors are indebted to the IDR/UPM Institute Angel Sanz-Andres for its kind support regarding the research related to power systems at the research institute. The authors are also grateful to Fernando Gallardo and SIELEC for the kind supervision and advise in relation to power distribution networks and protection systems. Finally, the authors are grateful to the reviewers, whose comments helped in improving the present work.

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

#### **References**


## *Article* **Detection Criterion for Progressive Faults in Photovoltaic Modules Based on Differential Voltage Measurements**

**Luis Diego Murillo-Soto 1, \* ,† and Carlos Meza 1,2,†**


**Abstract:** PV modules may experience degradation conditions that affect their power efficiency and affect the rest of the PV array. Based on the literature review, this paper links the parameter variation on a PV module with the six most common degradation faults, namely, series resistance degradation, optical homogeneous degradation, optical heterogeneous degradation, potential induced degradation, micro-cracks, and light-induced degradation. A Monte Carlo-based numerical simulation was used to study the effect of the faults mentioned above in the voltage of the modules in a PV array with one faulty module. A simple expression to identify faults was derived based on the obtained results. The simplicity of this expression allows integrating the fault detection technique in low-cost electronic circuits embedded in a PV module, optimizer, or microinverter.

**Keywords:** fault location in photovoltaic arrays; failure modes simulation; fault detection criterion

**Citation:** Murillo-Soto, L.D.; Meza, C. Detection Criterion for Progressive Faults in Photovoltaic Modules Based on Differential Voltage Measurements. *Appl. Sci.* **2022**, *12*, 2565. https://doi.org/10.3390/ app12052565

Academic Editors: Luis Hernández-Callejo, Maria del Carmen Alonso García and Sara Gallardo Saavedra

Received: 1 February 2022 Accepted: 18 February 2022 Published: 1 March 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/).

#### **1. Introduction**

According to [1], a fault can be defined as "an abnormal condition that may cause a reduction in, or loss of, the capability of a functional unit to perform a required function". Faults in photovoltaic (PV) installations can significantly affect their energy yield, which is why the development of fault detection and diagnosis techniques has become an essential topic of research in recent years, Refs. [2–4]. Research works that analyze failure mode in solar modules, e.g., Refs. [5–11], have identified that the power reduction in PV modules under fault causes a deformation of the current-voltage curve, e.g., curves in [12].

Recently, Ref. [13] presented a performance analysis with 30 faulty modules of different brands. The study shows that, on average, the power loss after two years, with several faults, drops around 1.08%; also, the voltage variation at maximum power was, on average −1.17%. Similar research in India concluded that, on average, power degradation in PV modules is found to increase 1.4% per year over 25 years [14]. But this degradation rate is not a static value; it depends on the life stage of the PV installations as is shown in [15], the first 18 years around 0.1% per year; then the next period of 10 years, the rate raises around 1%/year, and the final life stage, the rate increases more than 1.2%/year.

It is important to highlight that failure modes at the module level affect the whole PV installation performance because the PV array is a symmetrical composition of PV modules as is demonstrated by Gokmen in [16]. Eighteenth types of faults in solar modules have been described in [17], their detectable effects sometimes overlap, so it is not easy to distinguish among them. Furthermore, the diagnoses of faults require combining knowledge of different domains (visual inspection, thermography, electrical, chemical, material analyses, so on), Refs. [18,19]. This means that the specific failure modes could have, for example, the same behavior in the electric domain.

The detection criterion used in this paper is based on the widely accepted five parameters model, which, according to [20], can represent a single solar cell or several cells connected in series (PV module). This model is mathematically represented as follows,

$$I = I\_{ph} - I\_s \left( e^{\frac{V + IR\_s}{N\_s \eta V\_t}} - 1 \right) - \frac{(V + IR\_s)}{R\_p} \tag{1}$$

where *I* and *V* represent the current and voltage of the module, *Iph* is the photo-current generated, *I<sup>s</sup>* is the saturation current of the diode, *R<sup>s</sup>* and *R<sup>p</sup>* are the parasitic series and parallel resistances, *N<sup>s</sup>* is the number of cells in series, *η* represent the ideality factor, and *Vt* is the thermal voltage.

An increase in the PV module series resistance can model several fault mechanisms, as it has been indicated in [21,22]. However, using other parameters to explain other failure modes is scarce in the literature. In this regard, one of the contributions of this work is to propose parameter ranges in the PV model that represent others failure modes not included in the typical models, e.g., only series resistance variation. This work will focus on the most common types of degrading faults listed next.


Also, another PID fault variation, we call *PID*<sup>2</sup> , is found in [28], and they show that the parameters *R<sup>p</sup>* moves between [0.1, 100]%, *I<sup>s</sup>* varies between [100, 300]% and *η* moves between [100, 112]%. These percentages were calculated by analyzing the maximum and minimum values reported in the graphical results.


can produce changes of open-circuit voltage *Voc* in the range of [0, 18]% [20]. These percentages were calculated by analyzing graphical results.

A second contribution is the proposal of new detection criteria based on the measurements of the module voltages. Due to its simplicity, this criterion is suitable for implementation in monitoring strategies such as the one presented in [32]. In fact, several approaches use voltage measurements to detect and diagnose faults in the array, for example [33–38], but they require other variables like the string current, time, or even the ambient variables, to generate alarms adequately. Also, in those papers, the nature of faults analyzed is different because they focus mainly on faults such as diode short-circuit, open-circuits, partial shadows, and degradation faults represented only as an increase of the serial resistance. Furthermore, a similar indicator based on voltage is presented in [39], which is called ∆*V*, but this indicator requires knowing precisely the output voltage of the array and the theoretical voltage at its maximum power. This indicator is quite different from our approach since our criterion is calculated from the voltages of the modules in the PV string.

This work presents and analyzes new detection criteria to demonstrate that is capable to detect all the progressive faults only with the module voltage variable in multi-crystal PV modules. This numerical simulation study is divided into the following main sections: Section 2 describes the theoretical framework used to derive the detection criterion. Section 3 links specific fault conditions with changes in the variation of some parameters of the PV array. This section also explains the simulated experiment, circuits and software used, and the calibration process. Sections 4 and 5 show the main results and the discussion. Finally, the main conclusions are highlighted in section 6.

#### **2. Theoretical Framework: Detection Criterion Fundamentals**

The explicit expression of (1) for the module voltage is obtained using Lambert-W transformation, if we have an equation with the form *y* = *xe<sup>x</sup>* , the term *x* could be obtained as *<sup>x</sup>* <sup>=</sup> <sup>W</sup>(*y*) for any *<sup>x</sup>* <sup>&</sup>gt; <sup>−</sup>1/*e*. Therefore, the voltage is expressed as follows,

$$V = (I\_{ph} + I\_s)R\_p - I(R\_p + R\_s) - N\_s \eta V\_t \mathcal{W} \left(\frac{R\_p I\_s}{N\_s \eta V\_t} e \left(\frac{(I\_{ph} + I\_s - I)R\_p}{N\_s \eta V\_t}\right)\right). \tag{2}$$

It is clear that if a failure affects the PV module performance, the parameters in (1) and (2) change their values and hence the voltage *V* in (2) also varies [28]. If the initial values of all the parameters are known it is possible to identify the degradation condition of one PV module using (2). Nevertheless, the aforementioned is seldom the case. Not only the initial values are unknown but also a PV installation is comprised of several PV modules connected in series and parallel and the relationships between the parameter variation in one PV module will affect the others.

In a typical photovoltaic plant, PV modules are connected using a configuration known as series-parallel (SP). In an SP configuration, the PV modules are first connected in series forming strings which are then connected in parallel as depicted in Figure 1, forming a PV array of *m* × *n* modules. The voltages in the *j*-string, formed by *m* modules, are governed by Kirchhoff's voltage law [20], as follows;

$$\sum\_{i=1}^{m} V(i,j) - V\_{blk}(j) - V\_{op} = 0 \qquad , \tag{3}$$

where *V*(*i*, *j*) is the differential voltage of the module *i* in the string *j* as in Equation (2), *Vblk*(*j*) is the voltage of the blocking diode, and *Vop* is the operational voltage point of the array.

When all the modules are identical (3) can be simplified as follows,

$$mV(j) - V\_{blk}(j) - V\_{op} = 0.\tag{4}$$

Now consider the case in which at least one module is degraded or faulty. Due to the SP configuration the sum of all the voltages in the string with the faulty PV module does not change, see Equation (3). Therefore, the PV modules that are part of the string in which a faulty module is located will change their electrical operating point. We call the aforementioned elements *affected* PV modules, as it can be seen in Figure 2.

If we denote *m<sup>f</sup>* as the number of modules with the same fault, (3) can be rewritten as,

$$(m - m\_f) \cdot V\_a(j) + m\_f \cdot V\_f(j) = V\_{op} + V\_{blk}(j) \tag{5}$$

where *V<sup>f</sup>* (*j*) is the voltage of the faulty module (*degradeted*) and *Va*(*j*) is the voltage of the affected modules.

Notice that the sum of the voltages at the faulty and at the affected modules in a given string should be equal to the sum of the voltages of non-faulty modules in a string parallel to it. Hence, the variation in the voltages of the affected modules and the faulty modules must behave oppositely, i.e., when the voltage at the faulty decreases the voltage at the affected modules increases. Moreover, the difference between the voltages at the affected and the faulty modules can be expressed as follows,

$$
\Delta\_V = \frac{V\_a - V\_f}{V\_a} 100, \qquad \forall I > 0. \tag{6}
$$

Given that the current trough the string, *I*, is shared by the modules, (6) can be rewritten as

$$
\Delta\_V = \frac{IV\_a - IV\_f}{IV\_a} \\
100 = \frac{P\_a - P\_f}{P\_a} \\
100, \qquad \forall I > 0. \tag{7}
$$

where *P<sup>a</sup>* and *P<sup>f</sup>* are the power produced by the affected and faulty modules, and *I* is the current of the *j*-string (*j*). Therefore (6) represents a measure proportional to the power reduction in a PV module.

**Figure 1.** Nomeclature for the SP PV array.

**Figure 2.** Circuit simulated in LtSpice software with seven types of faults.

#### **3. Materials and Methods**

*3.1. Parameter Variation in Failure Modes*

The way in which photovoltaic modules degrade their power production is called failure mode, in other words, the process of how a PV module is going to fail is also known as a fault. The numerical simulation study considers variations in the parameters for 20 years according to the measured values of parameters reported in [17,20,23–30] for the faults mentioned in the Introduction. In this regard, Table 1 presents the range of variation of these parameters.


**Table 1.** Parameter variations for different fault families.

#### *3.2. Test Bench for Simulation*

The PV array was simulated using LtSpice [40] and was composed of 4 × 2 solar modules; each module uses the parameters of the Kyocera KC200GT PV module, as shown in Table 2. The simulated circuit is presented in Figure 2. The parameters used for the module are shown in Table 3 and were derived using the methodology proposed in [41]. Theoretical operative points are calculated and compared with the computing simulation to validate the general model of the 4 × 2 PV modules. The results are shown in Table 4 and this comparison is made with the relative error (*R*. *Error*) calculated as shown in (8).

#### *R*. *Error* = Theoretical value − Experimental value Theoretical value <sup>×</sup> 100. (8)

**Table 2.** Electrical Performance at Standard Test Conditions (STC).


**Table 3.** Values of the five-parameter model of the KC200GT [41].


**Table 4.** Electrical performance of the 4 × 2 array.


The simulations consist of selecting a PV module in the array and simulating each fault type 500 times. The selected PV module is called a faulty one, as shown in Figure 2. Furthermore, the numerical simulation changes the parameter values randomly according to Table 1 using a flat distribution over their predefined ranges. A sweep of voltage is made by the controlled voltage source for every new parameter setup, generating new curves. Also, it is important to mention that the voltage sweep varies from 0 to 4*Voc* V in steps of 0.5 V. In summary, every fault simulation has 500 parameter variations generating the following curves:


To simplify the experiments and for readability issues, the results for the seven simulated fault types are shown in Figure 3, and in the Appendix A from Figures A1–A6. These charts correspond only to faults located in the position (2, 2) in the array as indicated by the magenta rectangle in Figure 2. These curves can be extrapolated and are valid for faults in any position inside the array. All the simulations of the circuit were performed at standard test conditions, i.e., at an irradiance of 1000 W/m<sup>2</sup> and a cell temperature of 25 ◦ C, which are the main input references to get the electrical measurements reported in data-sheets [42]. Also, the datasheets give data to contrast with the simulation. The translation equations

were implemented according to [20], but no variations for these experiments were made on the temperature or the irradiance variables.

**Figure 3.** Series resistance degradation faults. This figure shows four charts as follows: (**a**) the current of the affected string versus the operating voltage of the array, (**b**) the differential voltage for all the modules in the affected string, (**c**) the relative percentage difference versus the operating voltage of the array, (**d**) peak behavior of output power in the array versus the operating voltage of the array.

#### **4. Results**

The results of the SRD fault are presented in Figure 3. The SRD fault was simulated by varying the resistance R4 up to ten times its original value. Figure 3a shows how the current of the faulty string changes as the degradation process caused by resistor R4 continues. This means that the differential voltages *V<sup>a</sup>* and *V<sup>f</sup>* are separating as the resistor R4 is increasing; see Figure 3b. The relative percentage difference calculated with Equation (6) is presented in Figure 3c. In this case, it is shown that there is a power loss by the faulty module concerning the non-faulty modules in the same string; the arrow shows the degradation process. The array's output power is shown in Figure 3d, and it goes down according to the increment of the serial resistance R4; also, the degradation process is indicated with the arrow.

Table 5 shows how the resistor value affects the voltage at the maximum powerpoint (*Vmpp*), the maximum power point (*Pmax*), and the total power loss (TPL). The TPL indicates how much power is lost in the whole PV array with respect to an expected non-faulty value of 1600 W, as indicated in Equation (8). Additionally, the table shows the value of voltage for the affected (*Va*) and faulty (*V<sup>f</sup>* ) PV modules as well as the relative percentage difference (∆*V*) between them. The last column in the table shows the module power loss (MPL) by the faulty module with respect to its maximum power as reported in the datasheet (200 W).

The graphical results for the following six fault types are presented in the Appendix A, each figure belongs to one fault type, and it contains the same four charts explained before for the SRD fault. In the Appendix A, the Figure A1 presents OD fault, and the results were obtained moving I4 values randomly. Figure A2 presents OHD fault; this simulation changes the values for R8 and I4 randomly and separately. Figure A3 presents PID<sup>1</sup> fault; here, the elements R8 and I4 move randomly but proportionally. Next, Figure A4 presents PID<sup>2</sup> fault; for this simulation, R8, I4, and the parameter *η* in diode D4 move freely and randomly. Further, the results for MC fault are presented in Figure A5, here I4, and the parameter *I<sup>s</sup>* in the diode D4 move randomly but in a proportional way. Finally, Figure A6 shows the results for LID fault when the parameter *I<sup>s</sup>* moves in the diode D4. It is interesting to see that all fault types present similar patterns and behaviors in all the charts.


**Table 5.** Numeric results for the increment of R4 in the 4 × 2 array.

#### **5. Discussion**

Only the PV module under examination has one type of progressive fault in the previous circuit simulations. The rest of the elements, such as PV modules, wires, bypass diodes, were assumed to work in a non-faulty condition. The aforementioned does not mean that the results are specific to the selected module location; if the faulty module is located at any other place in the array, the obtained charts would be equivalent to the graphs presented in this work.

Based on the simulation results it is confirmed that when a progressive fault affects a photovoltaic module, its differential voltage changes according to the severity of the failure mode. Figure 3b confirms that the progressive fault unbalances the voltage in the string; this means that as the progressive fault progress, the faulty module's differential voltage decrease, and the differential voltage of the affected modules increase. Hence, the relative percentage difference (∆*V*) increases, and the output power of the array decreases; these facts are easily checked in Figure 3c,d. For all the fault types analyzed, the same performance was observed in sub-charts in Figures A1–A6.

The degree of fault affectation in the module is calculated with the relative percentage difference (6), and this indicator is equivalent to the percentage of power loss if the modules work at the same string current. Table 5 presents the numerical results for the SRD fault. These results correspond to the evolution of the maximum power-point. In these maximum power points, it seems that the power loss indicator (∆*V*) is highly correlated with the power loss of the module calculated at standard test conditions, MPL. Actually, for this fault condition, the estimated power ∆*<sup>V</sup>* and the MPL are highly correlated, the Pearson Coefficient [43] is, in this case, *r* = 0.9999; therefore, it is possible to calculate a linear regression as,

$$MPL = 1.0311\Delta\_V - 0.1082,\ \forall \ \Delta\_V \text{ at } \text{Pmax.} \tag{9}$$

Moreover, the obtained linear regression has a very low variability, given a determination coefficient *R* <sup>2</sup> of 99.99%. This result is not unique for this failure mode, the correlation between ∆*<sup>V</sup>* and MPL appears for all the fault types analyzed. Graphically, it can be appreciated in Figure 4 and the parameters for the best fit line equation presented in Table 6.

In this PV array model of 1.6 kWp, as the progressive fault increase, the *Vmpp* moves affecting all the modules' power production; this means generating small losses in all the modules in the array. These small changes in power production caused by the progressive fault could be hidden because they are too small to be detected by the monitoring system. For the 1.6 kWp example, if *R<sup>s</sup>* double its value in whatever PV module, the total power loss of the array drops about 0.86%, which could be negligible. However, this percentage represents one faulty module losing around 6.26% of its energy (See the third row in Table 5). This hidden effect is even more drastic in large PV installations; let's suppose that a faulty module has a constant power loss, and if the number of PV modules in the array is high enough, the total power loss tends to zero. However, the faulty module is still losing the same percentage of energy and may evolve into a more severe fault condition.


**Table 6.** Parameters of the best fit line equation for all fault types.

The correlation *MPL* ∝ ∆*<sup>V</sup>* can be used to detect abnormal behavior. Moreover, it is possible to develop a simple criterion for fault detection and location in SP arrays. The Equation (10) is based on (6), and it checks if the relative percentage difference between two voltages is more than a threshold value called *δ*. Here *Vmax* represents the maximum differential voltage of all the modules in the string *j*, and *V*(*i*) is the differential voltage of the analyzed module. It should be noticed that *Vmax* is the higher voltage in the string, and it belongs to the module with less degradation. This criterion is simple to incorporate into other fault detection proposals that use differential voltages as input signals such as [44,45].

**Figure 4.** Power loss of the modules versus ∆*<sup>V</sup>* at mpp for the studied faults.

$$\text{Detect}(V\_{\text{max}}, V(i), \delta) = \begin{cases} \text{True} & \text{if } \quad \frac{V\_{\text{max}} - V(i)}{V\_{\text{max}}} > \delta \\\\ \text{False} & \text{if } \quad \frac{V\_{\text{max}} - V(i)}{V\_{\text{max}}} \le \delta \end{cases}, \forall i = \{1, 2, \dots, m\}. \tag{10}$$

To define the fault threshold, it is crucial to know the environment in which the solar array is placed. For instance, soiling may not be a problem in tropical regions benefiting from rainfall cleaning. On the contrary deserts or dry places, the PV array could be affected by dust storms or air pollution that could reduce the general performance [46,47]. In [48] it has been reported that power losses in outdoor conditions could be reduced between [5, 6]% and for laboratory conditions, it is possible to reduce it up to 40%. For instance, in tropical weather like Phitsanulok Tayland, it has been reported in [49] a decrease in solar radiation of [3.71, 11.15]% when the dust deposition rate (DDR) is 425 mg/m2d in 60 days. On the contrary dry cities like Mexico City for also 60 days, the DDR reported is 102 g/m2d to reduce the performance ratio up to 15%.

Detecting permanent faults such as PID, hotspots, or micro-cracks with online realtime methods is always a challenge. For instance, if a hotspot is considered mild, its temperature is just 10 ◦C higher than the other parts of the cells; however, if the hot spot is considered severe, it presents a temperature higher by approximately 18 ◦C. These facts mean that power production could be reduced between 4% and 10% [50]. A similar analysis can be done with micro-crack; for instance, in a PV module formed with 60 cells, if just one cell has an inactive area of 25%, the power loss in the whole module is about 10% [30]. Therefore, a rule to avoid false faults due to normal soiling or other temporal issues should use a *δ* value between 5 and 10 % as a threshold.

#### **6. Conclusions**

This work has analyzed several progressive faults presented in Table 1 and has been able to conclude the following:


• The power loss per module is proportional to the ∆*<sup>V</sup>* if the PV system works at maximum power-point.

The criteria proposed to detect permanent and progressive faults are based on estimating the power degradation in the faulty PV module. This is done with the differential voltages of the modules. This new detection criterion is suitable for real-time online analysis in PV arrays, which is part of the additional work, to experimentally demonstrate the simplicity of this fault detection criterion, running on a real-time system with several faulty modules.

**Author Contributions:** Conceptualization, L.D.M.-S.; methodology, L.D.M.-S.; software, L.D.M.-S.; validation, L.D.M.-S. and C.M.; formal analysis, L.D.M.-S.; investigation, L.D.M.-S. and C.M.; resources C.M.; data curation, L.D.M.-S.; writing—original draft preparation, L.D.M.-S.; writing review and editing, C.M.; visualization, L.D.M.-S.; supervision, C.M.; project administration, C.M.; funding acquisition, C.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the scholarship program of the Costa Rica Institute of Technology and the VIE project 5402-1360-4201.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Special thanks to Giovanni Spagnuolo from the University of Salerno, Italy, for his comments on this work.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A. Curves for the Simulated Failure Modes**

**Figure A1.** Optical homogeneous degradation faults. This figure shows four charts as following: (**a**) the current of the affected string versus the operation voltage of the array, (**b**) the differential voltage for all the modules in the affected string, (**c**) the relative percentage difference versus the operating voltage of the array, (**d**) peak behavior of output power in the array.

**Figure A2.** Optical heterogeneous degradation faults. This figure shows four charts as following: (**a**) the current of the affected string versus the operation voltage of the array, (**b**) the differential voltage for all the modules in the affected string, (**c**) the relative percentage difference versus the operating voltage of the array, (**d**) peak behavior of output power in the array.

**Figure A3.** Potential induced degradation faults type 1. This figure shows four charts as following: (**a**) the current of the affected string versus the operation voltage of the array, (**b**) the differential voltage for all the modules in the affected string, (**c**) the relative percentage difference versus the operating voltage of the array, (**d**) peak behavior of output power in the array.

**Figure A4.** Potential induced degradation faults type 2. This figure shows four charts as following: (**a**) the current of the affected string versus the operation voltage of the array, (**b**) the differential voltage for all the modules in the affected string, (**c**) the relative percentage difference versus the operating voltage of the array, (**d**) peak behavior of output power in the array.

**Figure A5.** Micro cracks faults. This figure shows four charts as following: (**a**) the current of the affected string versus the operation voltage of the array, (**b**) the differential voltage for all the modules in the affected string, (**c**) the relative percentage difference versus the operating voltage of the array, (**d**) peak behavior of output power in the array.

**Figure A6.** Light induced degradation faults. This figure shows four charts as follows: (**a**) the current of the affected string versus the operation voltage of the array, (**b**) the differential voltage for all the modules in the affected string, (**c**) the relative percentage difference versus the operating voltage of the array, (**d**) peak behavior of output power in the array.

#### **References**


## *Article* **Distribution Grid Stability—Influence of Inertia Moment of Synchronous Machines**

## **Tomáš Petrík, Milan Daneˇcek, Ivan Uhlíˇr, Vladislav Poulek and Martin Libra \***

Department of Physics, Faculty of Engineering, Czech University of Life Sciences Prague, Kamýcká 129, 16500 Prague, Czech Republic; tom.petr11@seznam.cz (T.P.); milan.danecek@gmail.com (M.D.); uhliri@tf.czu.cz (I.U.); poulek@tf.czu.cz (V.P.)

**\*** Correspondence: libra@tf.czu.cz

Received: 19 November 2020; Accepted: 14 December 2020; Published: 18 December 2020

**Abstract:** This paper shows the influence of grid frequency oscillations on synchronous machines coupled to masses with large moments of inertia and solves the maximum permissible value of a moment of inertia on the shaft of a synchronous machine in respect to the oscillation of grid frequency. Grid frequency variation causes a load angle to swing on the synchronous machines connected to the grid. This effect is particularly significant in microgrids. This article does not consider the effects of other components of the system, such as the effects of frequency, voltage, and power regulators.

**Keywords:** angle swinging; grid frequency oscillations; electromechanical system; inertial masses; microgrids

#### **1. Introduction**

The idea of grid voltage having a coherent sinusoidal course with constant frequency does not correspond with the reality of phenomena occurring in the distribution grid. Accurate measurements disclose oscillations of grid frequency and an instantaneous phase of distribution grid voltage. Grid frequency oscillation has a random noise trend because it is only the response of an energy system to a casually changing daily demand energy diagram [1]. This effect is particularly significant in microgrids [2–4]. The issue of network stability has been studied in a number of other works (e.g., [5]).

The phenomenon of grid oscillations may be described in terms of statistical dynamics using the power spectral density *S*ω(Ω) of a distribution grid's angular frequency fluctuation. An example of this function characterizing the frequency spectrum of grid frequency deviations is shown in Figure 1.

A spectrum of grid frequency deviations has two components: (i) A slow component for Ω < 2s −1 (mean ∆*F* < 0.32 Hz), which is determined by the properties of rotational regulation in power plants as well as the properties of superior central frequency regulation. The interconnection of originally regional districts of the distribution grid into larger complex grids and, later, into the European central grid caused a reduction in the amplitude of very slow fluctuations. Then, this phenomenon continued, leading to the current state of zero steady-state error in the distribution grid frequency. (ii) A fast component for Ω > 2s −1 (mean ∆*F* > 0.32 Hz), which is determined by an angular elasticity of the supply system in the case of power changes. A fast component originates from the magnitude and phase changes of instant voltage as it decreases on a power line, and the resulting fast load changes cause transformer impedances in the distribution grid. The angular elasticity refers to the angular elasticity of the phasors of voltage.

Figure 1 shows some samples of the spectrograms of fluctuation in the angular frequency of the distribution grid. Modern distribution grids with digital control of frequency have zero difference on average from central frequency, but there are visible fluctuations of frequency and phasor angle in an area up to Ω > 2 s −1 (mean ∆*F* > 0.32 Hz).

**Figure 1.** Examples of power spectral density *S*ω(Ω) at frequency fluctuation Ω.

If a synchronous machine driving inertial masses is connected to the grid, the mechanical movement of a machine shaft must follow frequency and phase changes of a voltage phasor of the grid. A shaft follows the instantaneous grid voltage phase with an angle deviation, producing a load angle β in a synchronous machine. The instantaneous power consumed or supplied by a machine to a grid depends on the magnitude and orientation of the load angle β. The greater the inertial masses on the shaft are in comparison to the size of a synchronous machine, the greater the dynamic deviation of the following fast grid phase changes and load angle oscillations will be [6].

Since an electromechanical system (a synchronous machine and its inertial masses) has very low oscillation dumps, the load angles, excited by distribution grid frequency fluctuations, have a harmonic swinging character. In addition, they are accompanied by an undesirable overflow of energy between the grid and the inertial masses on the shaft.

#### **2. Description of a Modeled System**

The electromechanical swinging system consists of the rigidity of a synchronous machine magnetic field, the damping effects, and the moment of inertia, which are connected with a shaft as presented in Figure 2. The system in Figure 2 is described by the differential equation

$$
\dot{I}\ddot{\beta} + \dot{B}\dot{\beta} + k\beta = \frac{I}{P\_d}\dot{\omega},
\tag{1}
$$

where β is the load angle of a synchronous machine, ω is the grid frequency (s−<sup>1</sup> ), *I* is the overall moment of inertia connected to a shaft (kg·m<sup>2</sup> ), *B* is the torsion damping constant (N·m·s), *K* is the torsion rigidity of a synchronous machine (N·m), and *P<sup>d</sup>* is the number of pole pairs in a synchronous machine.

**Figure 2.** Electromechanical vibrating set of a synchronous machine.

For simplification, the following will be assumed:


The precision of the system is described by Equation (1) in Laplace's transformation [7,8]:

$$G(p) = \frac{pI}{P\_d \left(p^2 I + p\mathcal{B} + k\right)}.\tag{2}$$

The power spectral density of the load angle fluctuation is

$$\mathcal{S}\_{\beta}(\Omega) = \lim\_{n \to \infty} (G(p))^2 \mathcal{S}\_{\omega}(\Omega),\tag{3}$$

where *S*β(Ω) is the power spectral density of the load angle fluctuation and *S*ω*(*Ω*)* is the power spectral density of fluctuation for angular frequency fluctuations (s−<sup>2</sup> ).

$$S\_{\beta}(\Omega) = \frac{\Omega^2 I^2 S\_{\omega}(\Omega)}{P\_d \left[ \left(k - \Omega^2 I\right)^2 + \Omega^2 \mathcal{B}^2 \right]}. \tag{4}$$

The natural frequency of the system is Ω0:

$$
\Omega\_0 = \sqrt{\frac{\mathbf{k}}{\mathbf{I}}}.\tag{5}
$$

The coefficient of a relative system damping

$$a = \frac{B}{2\sqrt{kl}}\tag{6}$$

is very low for most synchronous machines; it is usually < 0.1, in which case the resonance is very selective. The system transmits just the natural frequency. A system response has almost a sinusoidal course of load density swinging with an amplitude β*A*. From Equation (4), the results for Ω = Ω<sup>0</sup> can be written as

$$S\_{\beta}(\Omega\_0) = \beta\_A^2 = \frac{I^2}{P\_d^2 B^2} = S\_{\omega}(\Omega\_0). \tag{7}$$

#### **3. Maximum Admissible Grid Frequency Oscillations**

If the response of a load angle's swing has to have an amplitude *A* as the maximum, the power spectral density *S*ω(Ω0) of the grid angle frequency must not be greater than *S*ω,*max*(Ω0)

$$\mathcal{S}\_{\omega,\max}(\Omega\_0) = \frac{P\_d^2 \, B^2}{I^2} \, \beta\_{A'}^2 \tag{8}$$

where Sω*,max(*Ω0*)* is the upper limit of the power spectral density of a grid angle frequency (s−<sup>2</sup> ). Substituting the variable Ω<sup>0</sup> from Equation (5) into Equation (8) for the parameter *I* yields

*S*ω,*max*(Ω0) = *P* 2 *d B* 2 *k* 2 β 2 *<sup>A</sup>* <sup>Ω</sup><sup>4</sup> 0 . (9)

The maximum size of a flying wheel for which the amplitude of the load angle's swing does not exceed the given value β*<sup>A</sup>* is thus determined.

Courses for *S*ω*,max*(Ω0) according to Equation (9) for the chosen β*<sup>A</sup>* (e.g., 0.1, 0.15, and 0.2 rad) are drawn in Figure 3. These courses constitute the upper limits of power spectral density in the grid angle frequency fluctuations *S*ω,*max*(Ω0) for the permitted β*A*. The courses of the upper limits are fitted with a scaled *I* according to Equation (5):

$$I = \frac{k}{\Omega\_0^2} \tag{10}$$

**Figure 3.** Influence of the fluctuation spectrum of the grid and inertial momentum *I* on the amplitude of angular oscillation.

The maximum value of the moment of inertia being placed on the shaft of a synchronous machine is given by the first intersection point (from the right) of the function *S*ω,*max*(Ω0) for the permissible maximum of the load angle oscillation amplitude β*<sup>A</sup>* with the course *S*ω(Ω) measured in the grid section at the considered time.

#### **4. Practical Verifying**

The synchronous machine MEZ-A 225 MO 4 with the parameters *P<sup>s</sup>* = 50 kVA, *n* = 1500 min−<sup>1</sup> , *k* = 1140 N·m, and *B* = 14 N·m·s was chosen as an example. MEZ-A 225 MO 4 is a synchronous hydroalternator with four poles, without a dumper, with a brushless exciter powered by an external electronic source with a constant current of 1.2 A. This corresponds to excitation at a nominal voltage of 3 × 400 V in the idle state.

This excitation current of the brushless exciter was kept constant throughout the experiment.

The shaft of the machine was connected by a rigid coupling with a steel flywheel with a diameter of 1300 mm mounted in its own bearings.

The moment of inertia *<sup>I</sup>* <sup>=</sup> 160 kg·m<sup>2</sup> was determined by calculation from the geometric dimensions and the catalog data. The damping constant *B* = 14 N·m·s was calculated from the attenuation of the transient during phasing into the distribution network.

The unit was started at synchronous velocity by using a friction coupling from a diesel engine.

During the experiment, the synchronous machine ran as a motor on the distribution network without any additional load. It was loaded by mechanical losses in the bearings and the mechanical power consumption of its own ventilation and exciter. The total load was estimated at 3 kW, which is less than 8% of the nominal power of the synchronous machine. This state was close to the idling state, which is theoretically discussed above.

Under the permissible load angle oscillation amplitude β*<sup>A</sup>* = 0.1 rad (equal to the overflow of active power *<sup>P</sup><sup>a</sup>* <sup>=</sup> 20 kW), the intersection point P presented the maximum value *<sup>I</sup>* <sup>=</sup> 160 kg·m<sup>2</sup> .

Figure 4 gives evidence of conformity to the calculated value of an average amplitude (*<sup>I</sup>* <sup>=</sup> 160 kg·m<sup>2</sup> ) with a recording of grid frequency fluctuation and the corresponding response of the load angle. This recording was measured from the machine during the experiment.

**Figure 4.** Random changes in grid frequency (top) and changes influenced by the oscillation of the load angle of a synchronous machine (bottom).

#### **5. Conclusions**

Grid frequency variation causes a load angle to swing in the synchronous machines connected to the grid. The grid frequency variation spectrum shape grows the inertial moment on the shaft of the synchronous machine while also growing the swinging amplitude. This is visible from the shown equations. The maximum inertia moment value was thus determined for safe swinging amplitude. The evidence of conformity to the calculated value of the moment of inertia is shown in Figure 4 as a

recording of grid frequency and the corresponding response of the load angle. The average value of the load angle is present in the required interval. The instantaneous value of the load angle is mainly present in the required interval; however, the average value is more significant.

The load angle in this safe interval enables the machine to be driven in a safe and reliable manner by external controllers. These controllers are not considered in the article. This limitation should be respected both in drives with tens and hundreds of kW as well as in drives with small synchronous machines and stepping motors, where large inertial masses have to be driven.

As was written before, this paper did not consider external controllers. However, they are often a major part of frequency grid stability issues. Frequency grid fluctuations cause load angle fluctuations, which in turn cause frequency grid fluctuations. This interaction is a closed circle. In closing, it is necessary to mention trends and possibilities in frequency stabilization. A few studies last year discussed renewable sources of energy. These renewable sources can reduce the quality of grid services [9]. However, with the appropriate combination, along with batteries, they can be used to stabilize the properties of the grid [9]. A small battery source connected to photovoltaic panels was described in [10]. These battery sources, or even larger ones, would therefore be appropriate to include and use in the grid. This is especially true for intelligent buildings that combine many energy sources such as synchronous generators and photovoltaic panels; these batteries can be used as more than just backup sources.

**Author Contributions:** Conceptualization, T.P., I.U.; investigation, T.P., I.U., M.L.; methodology, M.D., V.P.; supervision, M.L.; validation, V.P., M.D.; visualization, V.P., M.L.; writing—original draft, I.U., T.P.; writing—review and editing, T.P., M.L. All authors have read and agreed to the published version of the manuscript.

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

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

#### **Nomenclature**


#### **References**


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

© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Analytical Modeling of Current-Voltage Photovoltaic Performance: An Easy Approach to Solar Panel Behavior**

**José Miguel Álvarez 1 , Daniel Alfonso-Corcuera 1,2 , Elena Roibás-Millán 1,2 , Javier Cubas 1,2 , Juan Cubero-Estalrrich 1 , Alejandro Gonzalez-Estrada 1 , Rocío Jado-Puente 1 , Marlon Sanabria-Pinzón <sup>1</sup> and Santiago Pindado 1,2, \***


**Abstract:** In this paper, we propose very simple analytical methodologies for modeling the behavior of photovoltaic (solar cells/panels) using a one-diode/two-resistor (1-D/2-R) equivalent circuit. A value of *a* = 1 for the ideality factor is shown to be very reasonable for the different photovoltaic technologies studied here. The solutions to the analytical equations of this model are simplified using easy mathematical expressions defined for the Lambert W-function. The definition of these mathematical expressions was based on a large dataset related to solar cells and panels obtained from the available academic literature. These simplified approaches were successfully used to extract the parameters from explicit methods for analyzing the behavior of solar cells/panels, where the exact solutions depend on the Lambert W-function. Finally, a case study was carried out that consisted of fitting the aforementioned models to the behavior (that is, the *I-V* curve) of two solar panels from the UPMSat-1 satellite. The results show a fairly high level of accuracy for the proposed methodologies.

**Keywords:** solar cell; solar panel; parameter extraction; analytical; Lambert W-function; spacecraft solar panels; *I-V* curve; modeling

#### **1. Introduction**

Based on the installation of power plants over the last few decades, it can be seen that photovoltaic energy has emerged as a very important factor in policies relating to renewable energy sources [1–7]. As an energy source, solar panels have a significant competitive edge over other renewable energy sources due to their potential for dual-use at either the industrial or the domestic scale. Both possibilities have provided an impetus for the increasing growth in demand for photovoltaic systems [5,8–11].

Photovoltaic energy has also been demonstrated to be crucial for stand-alone power systems such as glow-in-the-dark lettering and illuminated signs on highways. One of the first such stand-alone industrial applications of photovoltaic technologies was in power generation in spacecraft. According to Rauschenbach, "The first solar cell array that successfully operated in space was launched on 17 March 1958, on board Vanguard I, the second U.S. earth satellite" [12].

Modeling the performance of power sub-systems for space missions is essential at the predesign stage. At these early configuration stages of a space mission, simple, fast simulations are required that are as accurate as possible. This trend has also been amplified by the increasing importance of Concurrent Design (CD) in industrial processes. According to the European Cooperation for Space Standardization (ECSS), CD is especially important

**Citation:** Álvarez, J.M.; Alfonso-Corcuera, D.; Roibás-Millán, E.; Cubas, J.; Cubero-Estalrrich, J.; Gonzalez-Estrada, A.; Jado-Puente, R.; Sanabria-Pinzón, M.; Pindado, S. Analytical Modeling of Current-Voltage Photovoltaic Performance: An Easy Approach to Solar Panel Behavior. *Appl. Sci.* **2021**, *11*, 4250. https://doi.org/10.3390/app11094250

Academic Editor: Luis Hernández-Callejo

Received: 13 March 2021 Accepted: 3 May 2021 Published: 7 May 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/).

in the early predesign phase of space missions. It is based on the parallelization of processes associated with different sub-systems and is carried out in a special working environment organized into Concurrent Design Facilities (CDFs), where design parameters are shared, and the information flow is organized thus that stable solutions can be reached as quickly as possible [13–17].

Research carried out at the IDR/UPM Institute on solar panels has been driven by the need for simple procedures to calculate the performance of these panels. This need has arisen in the context of mission predesign at the CDF [18,19], and particularly in relation to sub-system coupling effects, such as those that can be found when analyzing thermal control and power distribution in spacecraft.

The simulation of photovoltaic devices (solar cells/panels) is normally carried out by means of equivalent circuit models, which are defined by implicit mathematical expressions. However, problems arise in terms of:


There are multiple ways to fit the equation for the photovoltaic equivalent circuit models to experimental data [20–29]. All of these approaches can essentially be grouped into two different types: Numerical (in which the equation is fitted to a large dataset that represents the *I-V* curve of the photovoltaic device) and analytical (in which the parameters of the equation are extracted based on three points of the *I-V* curve: short circuit, maximum power, and open voltage).

The second issue described above, i.e., the calculation of the output current as a function of the output voltage once the parameters of the equation have been extracted, can only be solved by numerical approaches [30–46] (which include iterative solutions [47–49]), or by using of the Lambert W-function [50–55] (whose exact calculation requires a numerical approach).

Possible solutions to the aforementioned problems include explicit methods of photovoltaic modeling, which are based on explicit equations with parameters defined based on the characteristic points of the *I-V* curve. Some of these methods also require the use of the Lambert W-function to define the parameters [56].

In the present work, we explore two important aspects of the analytical approach to the one diode/two resistor (1-D/2-R) equivalent circuit model for solar panels:

	- # Obtain the value of the first parameter of the 1-D/2-R equivalent circuit model, i.e., the value of the resistance of the series-connected resistor, and;
	- # Solve the implicit equation of this model to derive the output current in relation to the output voltage.

Our simplified approach to the Lambert W-function is also proven to be a relatively powerful tool for extracting the parameters of some of the most accurate explicit equations for photovoltaic modeling that can be found in the literature: the models of Kalmarlar and Haneefa [58,59] and Das [60].

The aim of this paper is to define very simple methodologies to measure the performance of photovoltaic devices (solar cells/panels), thus that these can be implemented as part of more complex simulations. One example of this type of simulation is the coupled thermo-electrical modeling of space systems with ESATAN© [61]. It should also be emphasized that these methodologies may be useful tools for professionals within small

and medium-sized enterprises in the solar energy sector to allow them to estimate the performance of the systems they design. Additionally, the increasing use of photovoltaic generation in small grids [62,63] could increase the worth of the easy approximations to solar panels performance included in this paper.

This paper is organized as follows. A simple model of solar panels is described in Section 2, and simplified equations for solving the 1-D/2-R equivalent circuit model and the Lambert W-function are given in Section 3. A case study is also described and solved in Section 3, using the procedures described in the preceding sub-sections. Finally, conclusions are presented in Section 4.

#### **2. Modeling of Photovoltaic Systems**

#### *2.1. The 1-D/2-R Equivalent Circuit Model*

The most widely used method of modeling the performance of a solar cell/panel (based on its *I-V* curve, where *I* is the output current and *V* the output voltage) is an equivalent circuit based on one current source, one diode, and two resistors (series and shunt resistors), as shown in Figure 1. The equation that describes the performance of this model is as follows [57,64–66]:

$$I = I\_{pv} - I\_0 \left[ \exp\left(\frac{V + IR\_{\rm s}}{naV\_T}\right) - 1\right] - \frac{V + IR\_{\rm s}}{R\_{\rm sh}} \tag{1}$$

where *Ipv* is the photocurrent, *I*<sup>0</sup> is the reverse saturation current of the diode, *Rsh* is the resistance of the shunt resistor, *R<sup>s</sup>* is the resistance of the series-connected resistor, and *V<sup>T</sup>* is the thermal voltage:

$$V\_T = \frac{\kappa T}{q} \tag{2}$$

**Figure 1.** (**a**) *I-V* curves for three different solar cells; (**b**) 1-D/2-R equivalent circuit model for analyzing the behavior of a solar cell/panel.

In the above equation, *κ* = 1.38064852·10 <sup>−</sup><sup>23</sup> m<sup>2</sup> kg s <sup>−</sup><sup>2</sup> K−<sup>1</sup> is the Boltzmann constant, *T* is the temperature expressed in K, and *q* = 1.60217662·10 <sup>−</sup><sup>19</sup> C is the electron charge. In Equation (1), *a* is the ideality factor of the diode (in principle, the effect of temperature on this parameter can be left aside [64]), and *n* is the number of series-connected cells in the solar panel (obviously, *n* = 1 when studying the performance of a single cell).

The first problem of modeling a photovoltaic device with the 1-D/2-R equivalent circuit model lies in extracting the five parameters of the model: *Ipv*, *I*0, *a*, *R<sup>s</sup>* , and *Rsh*. Depending on the available information, which might be either:


It is possible to extract these five parameters either by numerically fitting Equation (1) to the *I-V* curve, or by using the three characteristic points, which in practice represent four conditions [67]: the first three of which force the equation to match the points (0, *Isc*), (*Vmp*, *Imp*) and (*Voc*, 0), and the fourth is related to the maximum power condition:

$$-\frac{\partial I}{\partial V}\bigg|\_{\left(V\_{mp}, I\_{mp}\right)} = \frac{I\_{mp}}{V\_{mp}}.\tag{3}$$

Modeling the performance of photovoltaic systems with the equivalent circuit described here can be very advantageous, as these circuits somehow preserve the physical processes of a pair formed of a current source and a p-n junction [68]. The effect of the series-connected resistor, *R<sup>s</sup>* , is mainly associated with power losses in the solder bonds and interconnections between cells, whereas the effect of the shunt resistor, *Rsh*, is associated with current leakages across the p-n junction [69,70]. In recent years, several reviews of the different procedures and techniques for the parameter extraction problem related to the 1-D/2-R model have been published [20,21,23,24,67,71–73]. The present work focuses on an analytical approach based on the use of information from the characteristic points to calculate the five parameters in Equation (1). The following equations can be derived [57]:

$$\frac{n a V\_T V\_{\rm imp} \left(2I\_{\rm mp} - I\_{\rm sc}\right)}{\left(V\_{\rm mp} I\_{\rm sc} + V\_{\rm oc} \left(I\_{\rm mp} - I\_{\rm sc}\right)\right) \left(V\_{\rm mp} - I\_{\rm mp} R\_{\rm s}\right) - n a V\_T \left(V\_{\rm mp} I\_{\rm sc} - V\_{\rm oc} I\_{\rm mp}\right)} = \exp\left(\frac{V\_{\rm mp} + I\_{\rm mp} R\_{\rm s} - V\_{\rm oc}}{n a V\_T}\right),\tag{4}$$

$$R\_{\rm sh} = \frac{\left(V\_{mp} - I\_{mp}R\_{\rm s}\right)\left(V\_{mp} - R\_{\rm s}\left(I\_{\rm sc} - I\_{mp}\right) - naV\_{\rm T}\right)}{\left(V\_{mp} - I\_{mp}R\_{\rm s}\right)\left(I\_{\rm sc} - I\_{mp}\right) - naV\_{\rm T}I\_{mp}}\tag{5}$$

$$I\_{pv} = \frac{R\_{sh} + R\_s}{R\_{sh}} I\_{sc} \tag{6}$$

$$I\_0 = \frac{(R\_{sh} + R\_s)I\_{sc} - V\_{oc}}{R\_{sh} \exp\left(\frac{V\_{oc}}{naV\_T}\right)}.\tag{7}$$

Three problems arise at this point:


Fortunately, the ideality factor can be estimated by taking a value within the range from *a* = 1 to *a* = 1.5 [74,75], and the value of the resistance of the series-connected resistor, *Rs* , can be derived from [64]:

$$R\_s = A(W\_{-1}(B\exp(\mathbb{C})) - (D + \mathbb{C})),\tag{8}$$

where *W*−<sup>1</sup> is the negative branch of the Lambert W-function (see Section 2.1). The variables *A*, *B*, *C,* and *D* are defined as:

$$A = \frac{naV\_T}{I\_{mp}} \,\prime \tag{9}$$

$$B = -\frac{V\_{mp}\left(2I\_{mp} - I\_{sc}\right)}{V\_{mp}I\_{sc} + V\_{oc}\left(I\_{mp} - I\_{sc}\right)},\tag{10}$$

$$\mathcal{C} = -\frac{2V\_{mp} - V\_{oc}}{n a V\_T} + \frac{V\_{mp} I\_{sc} - V\_{oc} I\_{mp}}{V\_{mp} I\_{sc} + V\_{oc} \left(I\_{mp} - I\_{sc}\right)} \tag{11}$$

$$D = \frac{V\_{mp} - V\_{oc}}{naV\_T}.\tag{12}$$

However, the problem of solving the implicit equation in (1) still remains after the extraction of parameters *a* and *R<sup>s</sup>* . Authors such as Peng et al. [50] have shown that this equation can also be solved using the Lambert W-function:

$$I = \frac{R\_{\rm sh} \left(I\_{\rm pv} + I\_0\right) - V}{R\_{\rm sh} + R\_{\rm s}} - \frac{n n V\_T}{R\_{\rm s}} \mathcal{W}\_0 \left(\frac{R\_{\rm sh} R\_{\rm s} I\_0}{n n V\_T (R\_{\rm sh} + R\_{\rm s})} \exp\left(\frac{R\_{\rm sh} R\_{\rm s} \left(I\_{\rm pv} + I\_0\right) + R\_{\rm sh} V}{n n V\_T (R\_{\rm sh} + R\_{\rm s})}\right)\right),\tag{13}$$

where *W*<sup>0</sup> is the positive branch of the Lambert W-function (see Section 2.1).

#### *2.2. Explicit Equations/Models as Alternative to the 1-D/2-R Equivalent Circuit Model*

The difficulties associated with solving the 1-D/2-R model implicit equation drove researchers to develop explicit equations for modeling the performance of solar cells/ panels [58,59,76–82]. These photovoltaic models give the output current as a function of the output voltage by using simple mathematical models or equations that can be easily solved without the need for advanced mathematical tools. Of these, three models appear to be fairly accurate in relation to the *I-V* curve. The model proposed by Kalmarkar and Haneefa [58,59] is:

$$\frac{I}{I\_{\rm sc}} = 1 - (1 - \gamma)\frac{V}{V\_{\rm oc}} - \gamma \left(\frac{V}{V\_{\rm oc}}\right)^m \tag{14}$$

where [56]:

$$\gamma = \frac{2\left(\frac{I\_{mp}}{I\_{\rm sc}}\right) - 1}{(m-1)\left(\frac{V\_{mp}}{V\_{\rm ac}}\right)^{m}}\,\,\,\tag{15}$$

$$m = \frac{W\_{-1} \left( -\left(\frac{V\_{\rm oc}}{V\_{mp}}\right)^{\frac{1}{K}} \left(\frac{1}{K}\right) \ln\left(\frac{V\_{mp}}{V\_{\rm oc}}\right) \right)}{\ln\left(\frac{V\_{mp}}{V\_{\rm oc}}\right)} + \frac{1}{K} + 1,\tag{16}$$

and:

$$K = \frac{1 - \left(\frac{I\_{mp}}{I\_{\rm sc}}\right) - \left(\frac{V\_{mp}}{V\_{\rm sc}}\right)}{2\left(\frac{I\_{mp}}{I\_{\rm sc}}\right) - 1} \tag{17}$$

the model proposed by Das [60] is:

$$\frac{I}{I\_{\rm sc}} = \frac{1 - \left(\frac{V}{V\_{\rmoc}}\right)^k}{1 + h\left(\frac{V}{V\_{\rmoc}}\right)}\tag{18}$$

where [56]:

$$k = \frac{\mathcal{W}\_{-1}\left(\left(\frac{I\_{mp}}{I\_{\infty}}\right)\ln\left(\frac{V\_{mp}}{V\_{\infty}}\right)\right)}{\ln\left(\frac{V\_{mp}}{V\_{\infty}}\right)},\tag{19}$$

$$h = \left(\frac{V\_{oc}}{V\_{mp}}\right) \left(\frac{I\_{sc}}{I\_{mp}} - \frac{1}{k} - 1\right),\tag{20}$$

while the model of Pindado and Cubas [56,80] is:

$$I = \begin{cases} \begin{array}{l} I\_{\rm SC} \left( 1 - \left( 1 - \frac{I\_{mp}}{I\_{\rm sc}} \right) \left( \frac{V}{V\_{mp}} \right)^{\frac{I\_{mp}}{1sc - I\_{mp}}} \right) \vdots \, V \le V\_{mp} \\\ I\_{mp} \left( \frac{V\_{mp}}{V} \right) \left( 1 - \left( \frac{V - V\_{mp}}{V\_{oc} - V\_{mp}} \right)^{\Phi} \right) \vdots \, V \ge V\_{mp} \end{array} \tag{21}$$

where:

$$\phi = \left(\frac{I\_{\rm sc}}{I\_{mp}}\right) \left(\frac{I\_{\rm sc}}{I\_{\rm sc} - I\_{mp}}\right) \left(\frac{V\_{\rm oc} - V\_{mp}}{V\_{\rm oc}}\right). \tag{22}$$

In the graph in Figure 2, these models are fitted to the experimental *I-V* curve for an RTC solar cell, taken from the well-known work by Easwarakhanthan et al. [83]. The 1-D/2-R equivalent circuit model is also fitted to these data, and it can be seen that the accuracy of the explicit models is quite high, as stated in [56,81]. The problem with using some of these models again lies in the need to solve the Lambert W-function, which requires some mathematical ability.

**Figure 2.** *I-V* curve of the RTC solar cell [83], with the results from the 1-D/2-R equivalent circuit model, and the explicit models of Kalmarkar and Haneefa, Das, and Pindado and Cubas.

#### *2.3. The Lambert W-Function*

From the methodologies described above, it is clear that in order to work with the equations of the implicit 1-D/2-R model and some of the explicit models, some knowledge of the Lambert W-function is required. The Lambert W-function, *W*(*z*) (plotted in Figure 3), is defined as:

$$z = W(z) \exp(W(z)),\tag{23}$$

In the above equation, *z* is a complex number. If a real variable *x* is considered, the Lambert function is then defined within the range [−1/e, <sup>∞</sup>]. It should also be noted that this function gives a double value within the range [−1/e, 0]. Two different branches, called the positive and negative branches, are normally defined for this function, as follows:


This function is a useful tool for solving equations that involve exponentials since if *X* = *Y*exp(*Y*), then *Y = W(X)*. However, as stated above, a certain level of expertise is required to solve these [84,85]. Barry et al. developed an interesting explicit approach, although this may not be sufficiently direct to obtain a solution [86].

**Figure 3.** The Lambert W-function, with its two different branches.

#### **3. Results**

In this section, a method that facilitates the use of the 1-D/2-R equivalent circuit model is presented. It is organized into three sub-sections, as follows:


#### *3.1. On the Best Value for the Ideality Factor*

In a paper published in the 2nd International Conference on Renewable Energy Research and Applications (ICRERA 2013), our research group demonstrated a hyperbolic relationship between the non-dimensional RMSE, *ξ*, of the 1-D/2-R equivalent circuit model fitted to the *I-V* curve of a solar cell, and the value of the ideality factor, *a* [87]. The non-dimensional (or normalized) RMSE is defined as:

$$\mathcal{L} = \frac{\text{RMSE}}{I\_{\text{sc}}} = \frac{1}{I\_{\text{sc}}} \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left( I - I\_{\text{ref}} \right)^{2}} \,\text{}\tag{24}$$

where *N* is the number of point of the *I-V* curve, *I* is the current obtained from modeling the corresponding solar cell/panel, and *Iref* is the measured current associated with each one of those points. The aforementioned relationship between *ξ* and *a* is plotted in Figure 4. It can be seen that the curve has a hyperbolic shape, with a minimum value of the nondimensional RMSE at *a* = 1.18. This hyperbolic shape was also confirmed in a recent work by Elkholy and Abou El-Ela [88].

The normalized RMSE can be used to compare the accuracy of different models applied to different photovoltaic technologies. However, it should also be taken into account that variations in temperature and irradiance could somehow affect the accuracy of the model's accuracy. In the present work, the current at the short circuit from the reference *I-V* curve has been selected as it was in previous works such as [57]. This way to normalize the RMSE has been used by other authors [89]. However, it was not the only one, as the average current value from the measured dataset that represents the *I-V* curve [88,90], or the difference between the current at two points [91], has also been proposed as current values to normalize the RMSE.

**Figure 4.** Non-dimensional RMSE, *ξ*, for a 1-D/2-R equivalent circuit model fitted to the *I-V* curve of an Emcore ZTJ solar cell (see Figure 1) versus the ideality factor, *a*.

The value of the ideality factor, *a*, was studied with the aim of providing future researchers with an accurate but much simpler fitting of the 1-D/2-R equivalent circuit model. To carry out this analysis, the *I-V* curves from four solar cells and three solar panels were selected from previous works [56,81], as shown in Table 1. The parameters extracted by fitting Equation (1) to these *I-V* curves using Matlab © are given in Table 2. These fits were used as a reference for the values calculated analytically since it could be reasonably assumed that they were the most accurate ones for the *I-V* curves. It should also be remarked that the equations for the 1-D/2-R equivalent circuit model resulting from these fits may not exactly meet the requirements for the short circuit, MPP, and open-circuit conditions, as they represent the best fits not only to the three characteristic points but to a large number of additional points.


**Table 1.** Solar cell/panel *I-V* curves used in the present work (where *n* is the number of series-connected cells for each of these photovoltaic devices; *T* is the temperature; and the characteristic points are *Isc, Imp, Vmp*, and *Voc*).

<sup>1</sup> Taken from Easwarakhantan et al. [83]. <sup>2</sup> Graphically extracted from the manufacturer's datasheet. <sup>3</sup> Supplied by Azur Space. <sup>4</sup> Measured at CIEMAT (Spain).

**Table 2.** Parameters for the 1-D/2-E equivalent circuit models, numerically fitted to *I-V* curves based on data from RTC France, TNJ Spectrolab, ZTJ Emcore, and Azrur Space 3G30C solar cells, and Photowatt PWP 201, Kyocera KC200GT-2, and Selex Galileo SPVS X5 solar panels (where *ξ* is the non-dimensional RMSE for the fit).


The 1-D/2-R equivalent circuit model was also analytically fitted to the aforementioned *I-V* curves. Using specific values for the ideality factor, *a*, and the characteristic points of the curve, *Isc*, *Vmp*, *Imp* and *Voc*, the different parameters of the model were found using the following procedure:


Once all the parameters were known, the *I-V* curve was found by obtaining the values of the current, *I*, for each value of the voltage, *V*, using Equation (13). The ideality factor was initially varied from *a* = 0.5 to *a* = 1.8 for each photovoltaic device (although this range was extended for some devices), and the non-dimensional RMSE, *x*, was calculated in each case. The results are plotted in Figure 5. Table 3 shows the values of the parameter *a* at which the minimum value of *x* is reached, with the rest of the parameters for the 1-D/2-R equivalent circuit.

**Figure 5.** Non-dimensional RMSE, *ξ*, corresponding to the analytical approximations for extracting the 1-D/2-R equivalent circuit parameters in relation to the ideality factor, *a*, selected, for each one of the analyzed solar cells and panels.

**Table 3.** Parameters of the 1-D/2-R equivalent circuit models analytically fitted to *I-V* curves corresponding to data from RTC France, TNJ Spectrolab, ZTJ Emcore, and Azrur Space 3G30C solar cells, and Photowatt PWP 201, Kyocera KC200GT-2, and Selex Galileo SPVS X5 solar panels (where *ξ* is the non-dimensional RMSE for the fit).


Some important conclusions can be drawn based on these results, as follows:


of the shunt resistor shown in Table 3, which is obtained for the Photowatt PWP 201 solar panel.

• The best analytical fit gave remarkably low values of a in two cases: the Azur Space 3G30C solar cell, and the Selex Galileo SPVS X5 solar panel (composed of five Azur Space 3G28C series-connected cells), as their values were outside the range [1,1.5] suggested by Villalva et al. [74,75].

In order to compare both types of fit, the values of the ideality factor that gave the lowest values of x for the analytical fittings, *abf-a* , were plotted against the values obtained from the numerical fittings, *abf-n*, in Figure 6. It can be seen that the correlation between the fits is high, with the exception of only two devices: the Azur Space 3G30C solar cell, and the Selex Galileo SPVS X5 solar panel (which are both based on the same photovoltaic technology, as shown in Table 1). However, this correlation does not allow us to identify any rule that would suggest a reasonable value for the ideality factor. We, therefore, plotted the average values of the non-dimensional RMSE, *ξav*, obtained from the analytical fits versus the ideality factor, *a*. It can be seen from Figure 6 that the minimum value of *ξav* was obtained at *a* = 0.98. As a consequence, the value suggested for all devices is *a* = 1.

**Figure 6.** (**a**) Ideality factor corresponding to the best analytical fit to the 1-D/2-R equivalent circuit model, *abf-a* , for the studied *I-V* curves, versus ideality factors obtained from the numerical fittings, *abf-n* . The line corresponding to perfect correlation between both factors is shown as a reference. (**b**) Averaged non-dimensional RMSE for all analytical fittings, *ξav*, versus the ideality factor, *a*, used in the fits.

Figure 7 shows a comparison of the non-dimensional RMSE values obtained from the numerical fit, the analytical fit, and the analytical fit for *a* = 1 (see Table 4) for all of the photovoltaic devices. Although the differences may appear high, it should be taken into account that photovoltaic systems are normally designed to operate at the MPP, which is exactly reproduced by all of the analytical approximations for solving the parameter extraction of the 1-D/2-R equivalent circuit model. In addition, we reviewed the results from several works published over the last five years on parameter extraction for 1-D/2-R equivalent circuit models, based on the same experimental data from the RTC France solar cell and the Photowatt PWP 201 solar panel. The values obtained in these works for the non-dimensional RMSE for the fits were lower than the results from our method, with average values of *ξ* = 1.20·10 −03 (RTC France solar cell) and *ξ* = 4.03·10 −03 (Photowatt PWP 201 solar panel). Although these more accurate values were based on numerical rather than analytical procedures, it should be pointed out that numerical approaches do not always reach a more accurate solution [92].

**Figure 7.** Values for the non-dimensional RMSE, *x*, obtained from a numerical fit of the 1-D/2-R equivalent circuit model to the studied *I-V* curves (see Table 2), an analytical fit (see Table 3), and an analytical fit with *a* = 1 (see Table 4).

**Table 4.** Analytically fitted parameters for 1-D/2-R equivalent circuit models, with *a* = 1, for *I-V* curves corresponding to RTC France, TNJ Spectrolab, ZTJ Emcore, and Azrur Space 3G30C solar cells, and Photowatt PWP 201, Kyocera KC200GT-2, and Selex Galileo SPVS X5 solar panels (where *ξ* is the non-dimensional RMSE for the fit).


#### *3.2. Lambert W-Function Simplified Equations for Solar Cell/Panel Modeling*

A thorough review of the available literature from between 2000 and 2020 was carried out to find relevant data on photovoltaic devices, including the five parameters of the 1-D/2-R equivalent circuit model (*Ipv*, *I*0, *a*, *R<sup>s</sup>* , and *Rsh* ), the characteristic points of the *I-V* curve (*Isc*, *Vmp*, *Imp*, and *Voc*), the temperature of the cells, *T*, which is related to the performance curve, and the number of series-corrected cells, *n*. Information was found for 90 photovoltaic devices, most of which were solar panels.

The positive branch of the Lambert W-function needs to be solved for in Equation (13), and must be evaluated at certain points x for a given value of the output voltage, *V* (within the range [0, *Voc*]):

$$\mathbf{x} = f\left(I\_{\mathrm{pv}}, I\_{0\prime}a, \mathcal{R}\_{\mathrm{s}}, \mathcal{R}\_{\mathrm{sh}}, \mathbf{n}, T, V\right). \tag{25}$$

After post-processing the data from the 90 solar cells/panels, it was clear that Equation (13) refers to the right-side section of the Lambert's W-function positive branch, W<sup>0</sup> + . The following expressions were proposed:

$$\mathcal{W}\_0^+(\mathbf{x}) = \mathbf{x} \exp\left(0.71116 \mathbf{x}^2 - 0.98639 \mathbf{x}\right), \mathbf{x} \in \left[2 \cdot 10^{-16}, 2 \cdot 10^{-1}\right],\tag{26}$$

$$\mathcal{W}\_0^+(\mathbf{x}) = -1.6579 + 0.1396 \left( 2.9179 \cdot 10^5 - \left( \mathbf{x} - 22.8345 \right)^4 \right)^{0.25}, \mathbf{x} \in [0.2, 1.2], \tag{27}$$

$$\mathcal{W}\_0^+(\mathbf{x}) = -1.2216 + 3.4724 \cdot 10^{-2} \left( 1.7091 \cdot 10^8 - \left( \mathbf{x} - 114.146 \right)^4 \right)^{0.25}, \mathbf{x} \in [1.2, 10]. \tag{28}$$

The error for these equations was below 0.15% (Equation (26)), 0.16% (Equation (27)), and 0.08% (Equation (28)). For larger values of *x*, we recommend using the approximation proposed by Barry et al. [86] (see Appendix A). In Figure 8, the positive branch of Lambert's W-function, *W*<sup>0</sup> + , calculated at the points *x* (Equation (25)) corresponding to the data for the solar cells/panels, evaluated at *V* = 0 and *V* = *Voc*, is shown together with Equation (26) to (28). Depending on the magnitude of *x* (that is, depending on how closely it approaches 0 + ), *W*<sup>0</sup> + can be reasonably estimated as:

$$W\_0^+ = \mathbf{x} - 0.98639\mathbf{x}^2 + 1.1976\mathbf{x}^3 \approx \mathbf{x} - 0.98639\mathbf{x}^2 \approx \mathbf{x},\tag{29}$$

The negative branch of the Lambert W-function, *W*−1, needs to be calculated at:

$$\mathfrak{x} = f\left(I\_{p\upsilon\nu}I\_{0\prime}a, \mathcal{R}\_{\text{s}\prime}\mathcal{R}\_{\text{sh}\prime}\mathfrak{n}, T, V\right) \tag{50}$$

when estimating the series resistance, Rs, for the 1-D/2-R equivalent circuit model (Equation (8) to (12)). Since the above equation for the different solar cells/panels from the database gives values of *x* within the range [−2.77·10 −3 , −10 −35 ], two expressions were defined:

$$\begin{array}{l} W\_{-1}(\mathbf{x}) = 9.7117 \cdot 10^{-5} \ln(-\mathbf{x})^3 + 6.8669 \cdot 10^{-3} \ln(-\mathbf{x})^2 + 1.2 \ln(-\mathbf{x}) - 1.1102, \\ \mathbf{x} \in \left[ -10^{-2}, -5 \cdot 10^{-13} \right] \end{array} \tag{31}$$

$$\begin{cases} \mathcal{W}\_{-1}(\mathbf{x}) = 1.6705 \cdot 10^{-6} \ln(-\mathbf{x})^3 + 4.4514 \cdot 10^{-4} \ln(-\mathbf{x})^2 + 1.0511 \ln(-\mathbf{x}) - 2.3364 \\ \mathbf{x} \in \left[ -5 \cdot 10^{-13}, -10^{-40} \right] \end{cases} \text{. (32)}$$

The error for these equations was below 0.42% (Equation (31)) and 0.02% (Equation (32)) (see Figure 9). Equation (32) can be applied for values of up to *x* = −10 −50 , with errors of below 0.07%.

**Figure 9.** (**a**) Negative branch of Lambert's W-function, W-1 , calculated for points *x* (Equation (30)) corresponding to the series resistance, *Rs*, of the 1-D/2-R equivalent circuit model, based on data for solar cells/panels. Our approximate expression for *W*−<sup>1</sup> (Equations (31) and (32)) is also shown. (**b**) Negative branch of Lambert's W-function, *W*−<sup>1</sup> , calculated for points *x* (Equation (33)) corresponding to the explicit models proposed by Kalmarkar and Haneefa and Das, based on data for solar cells/panels. Our approximate expression for *W*−<sup>1</sup> (Equations (35) and (36)) is also shown.

Finally, if the explicit models proposed by Kalmarkar and Haneefa, and Das are used, W-1 needs to be calculated at points x defined by the characteristic points of the *I-V* curve (*Isc*, *Vmp*, *Imp*, and *Voc* (see Equations (16) and (19)):

$$\mathfrak{x} = f\left(I\_{\mathrm{sc}\prime}I\_{mp\prime}V\_{\mathrm{oc}\prime}V\_{mp}\right) \tag{33}$$

as shown in Figure 9. Based on the post-processed data for the 90 solar cells/panels, the values of *x* calculated using these explicit methods were within the range [−0.304, −0.1]. The following equation for *W*-1 was initially suggested, which was characterized by an error of less than 1.6%22 within this range:

$$\mathcal{W}\_{-1}(\mathbf{x}) = 248.42\mathbf{x}^4 + 134.24\mathbf{x}^3 + 4.4258\mathbf{x}^2 - 14.629\mathbf{x} - 4.9631\tag{34}$$

In order to extend both the range of validity and the accuracy of the above equation, a splitting into two sub-branches at the point where the curvature of the function changes its sign (*x* = −0.27), the following equations were suggested:

$$\mathcal{W}\_{-1}(\mathbf{x}) = -1 - \sqrt{42.949\mathbf{x}^2 + 37.694\mathbf{x} + 8.0542}, \ \mathbf{x} \in [-0.36785, -0.27] \tag{35}$$

$$\begin{array}{l} \mathcal{W}\_{-1}(\mathbf{x}) = 0.14279 \ln(-\mathbf{x})^3 + 1.04416 \ln(-\mathbf{x})^2 + 3.92 \ln(-\mathbf{x}) + 1.65795, \\ \mathbf{x} \in [-0.27, -0.0732] \end{array} \tag{36}$$

The error for the above equations was below 0.23% (Equation (35)) and 0.02% (Equation (36)) (see Figure 9). It should also be emphasized that caution needs to be used with regard to Equation (35), as it gives complex solutions when *x*→−1/*e* + (that is, when *x* ∈ [−1/*e*, −0.36785]).

#### *3.3. Case Study: The Solar Panels of the UPMSat-1*

The general characteristics from the *I-V* curves for the UPM-5 and UPM-6 solar panels for UPMSat-1 (see Figure 10) are shown in Tables 1, 2 and 5 (Appendix B). These panels are based on two different technologies: Si (UPM-5) and Ga-As (UPM-6). In the latter, the number of series-connected cells *n* was 32; however, since they were based on double-junction technology, this number must be multiplied by two, unlike in the

1-D/2-R equivalent circuit model (the Selex Galileo SPVS X5, formed from five Azur Space triple-junction technology series-connected solar cells, where *n* = 15, as shown in Table 1).

**Figure 10.** The UPMSat-1 satellite; launched in 1995, this was the tenth university-developed space mission in history [93]. (**a**) Satellite during the integration in the Ariane IV-40 launcher (V75 flight). (**b**) Exploded view drawing.

**Table 5.** Characteristics of the measured *I-V* curves corresponding to the UPM-5 and UPM-6 solar panels of UPMSat-1 (see Figure 10).


The *I-V* curves for the UPM-5 and UPM-6 solar panels are plotted in Figures 11 and 12, respectively, with several other curves corresponding to:


**Table 6.** Parameters for the 1-D/2-R equivalent circuit model (where Num. represents the numerical best fit, and An. represents the analytical fit with *a* = 1) for the UPM-5 and UPM-6 solar panels of the UPMSat-1 satellite. *ξ* is the non-dimensional RMSE for the experimental data.



**Table 7.** Parameters for the explicit models proposed by Kalmarkar and Haneefa, Das and Pindado and Cubas, for the UPM-5 and UPM-6 solar panels of the UPMSat-1 satellite. *ξ* is the non-dimensional RMSE for the experimental data.

**Figure 11.** *I-V* experimental curve (left axis) for the UPM-5 solar panel. Curves for the 1-D/2-R equivalent circuit model (obtained using two procedures: best fit to the parameters obtained numerically and analytical extraction) and the explicit methods proposed by Kalmarkar and Haneefa, Das and Pindado and Cubas are shown. The differences in the current compared to the experimental data are indicated on the right axis.

**Figure 12.** *I-V* experimental curve (left axis) for the UPM-6 solar panel. Curves for the 1-D/2-R equivalent circuit model (obtained using two procedures: best fit to the parameters obtained numerically and analytical extraction) and the explicit methods of Kalmarkar and Haneefa, Das and Pindado and Cubas are shown. The differences in the current compared to the experimental data are indicated on the right axis.

The differences in the current for these models with regard to the experimental *I-V* curves are plotted in Figures 11 and 12.

The differences were relatively small, with the maximum always located in a reasonably small range between the MPP and the open voltage point. The values of the non-dimensional RMSE, *ξ*, for the fits are plotted in Figure 13. The results show reasonably good agreement with the mathematical approaches proposed in the present work.

**Figure 13.** Non-dimensional RMSE, *ξ*, for the numerical fit of the 1-D/2-R equivalent circuit model, and the analytical fit with the proposed method, for the UPM-5 and UPM-6 solar panels of UPMSat-1. Values derived from the explicit methods studied here are also shown.

#### **4. Conclusions**

In this study, we have reviewed analytical approaches for extracting the five parameters of the 1-D/2-R equivalent circuit model from the characteristic points of the *I-V* curve (short circuit, MPP, and open circuit points). Five different problems were identified, as follows:


*I-V* curves from seven different solar cells/panels were used to test our proposed approach. *I-V* curves from two of the solar panels of the UPMSat-1 spacecraft were also analyzed as a case study.

The most important conclusions of this work are as follows:


Our approach was carefully verified in a case study in which the fits of the1-D/2-R equivalent circuit model and the explicit methods to the measured data (*I-V* curve) for two solar panels from the UPMSat-1 satellite were compared with the results of the proposed method.

Finally, it should be highlighted that the results from this work open up new possibilities for coupled calculations, both for the performance of photovoltaic systems and in other disciplines such as thermodynamics and power distribution in grids. The simple but accurate solutions to the 1-D/2-R equivalent circuit model and the explicit methods described here can easily be implemented in software packages such as ESATAN©.

**Author Contributions:** Project administration, S.P., J.C., E.R.-M.; methodology, S.P.; conceptualization, S.P.; formal analysis, S.P.; supervision S.P., J.C., E.R.-M.; investigation, J.M.Á., D.A.-C., J.C.-E., M.S.-P., R.J.-P., A.G.-E.; software, J.M.Á., D.A.-C., J.C.-E., M.S.-P., R.J.-P., A.G.-E.; validation, visualization, J.M.Á., D.A.-C., J.C.-E., M.S.-P., R.J.-P., A.G.-E.; writing—original draft, J.M.Á., D.A.-C., E.R.-M., J.C.-E., M.S.-P., R.J.-P., A.G.-E.; writing—review and editing, S.P. 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:** This work was carried out within the framework of the Educational Innovation Project PIE IE1920.1402 'PIRAMIDE: *Proyectos de Investigación Realizados por Alumnos de Máster/Grado para la Innovación y el Desarrollo Espacial*'. This Academic Innovation Project (AIP) was funded by *Universidad Politécnica de Madrid* (UPM) in the 2019/2020 call '*Ayudas a la innovación educativa y a la mejora en la calidad de la enseñanza*'. The authors are indebted to Ángel Sanz-Andrés for his support with regard to the research program on the performance of solar panels at *Instituto Universitario de Microgravedad "Ignacio Da Riva"* (IDR/UPM), *Universidad Politécnica de Madrid* (UPM). The authors are grateful to the reviewers for their kind help in improving the present paper.

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

#### **Appendix A. Approximation to the Right-Side of the Lambert W-Function Positive Branch by Larry et al. (2000) (see References)**

$$\mathcal{W}\_0^+(\mathbf{x}) = 1.4586887 \ln \left( \frac{1.2 \mathbf{x}}{\ln(2.4 \mathbf{x}/\ln(1 + 2.4 \mathbf{x}))} \right) - 0.4586887 \ln \left( \frac{2 \mathbf{x}}{\ln(1 + 2 \mathbf{x})} \right). \tag{A1}$$

#### **Appendix B.** *I-V* **Curves of the UPMSat-1 Solar Panels.**


**Table 1.** *I-V* curve of the UPM-5 solar panel of the UPMSat-1.


**Table 1.** *Cont.*

**Table 2.** *I-V* curve of the UPM-6 solar panel of the UPMSat-1.



**Table 2.** *Cont.*

#### **References**


**Andrés Pérez-González \* , Nelson Benítez-Montoya , Álvaro Jaramillo-Duque and Juan Bernardo Cano-Quintero**

> Research Group in Efficient Energy Management (GIMEL), Electrical Engineering Department, Universidad de Antioquia, Calle 67 No. 53-108, Medellín 050010, Colombia; nelson.benitez@udea.edu.co (N.B.-M.) alvaro.jaramillod@udea.edu.co (Á.J.-D.); bernardo.cano@udea.edu.co (J.B.C.-Q.) **\*** Correspondence: afernando.perez@udea.edu.co; Tel.: +57-322-639-7758

**Abstract:** Solar energy is one of the most strategic energy sources for the world's economic development. This has caused the number of solar photovoltaic plants to increase around the world; consequently, they are installed in places where their access and manual inspection are arduous and risky tasks. Recently, the inspection of photovoltaic plants has been conducted with the use of unmanned aerial vehicles (UAV). Although the inspection with UAVs can be completed with a drone operator, where the UAV flight path is purely manual or utilizes a previously generated flight path through a ground control station (GCS). However, the path generated in the GCS has many restrictions that the operator must supply. Due to these restrictions, we present a novel way to develop a flight path automatically with coverage path planning (CPP) methods. Using a DL server to segment the region of interest (RoI) within each of the predefined PV plant images, three CPP methods were also considered and their performances were assessed with metrics. The UAV energy consumption performance in each of the CPP methods was assessed using two different UAVs and standard metrics. Six experiments were performed by varying the CPP width, and the consumption metrics were recorded in each experiment. According to the results, the most effective and efficient methods are the exact cellular decomposition boustrophedon and grid-based wavefront coverage, depending on the CPP width and the area of the PV plant. Finally, a relationship was established between the size of the photovoltaic plant area and the best UAV to perform the inspection with the appropriate CPP width. This could be an important result for low-cost inspection with UAVs, without high-resolution cameras on the UAV board, and in small plants.

**Keywords:** deep learning (DL); unmanned aerial vehicle (UAV); photovoltaic (PV) plants; semantic segmentation; coverage path planning (CPP)

#### **1. Introduction**

According to REN21, over the past two years, global photovoltaic (PV) plants capacities and annual additions have grown and expanded rapidly. For instance, 621 GW were installed in the year 2019 and 760 GW in the year 2020 [1], despite the reduction in electricity consumption and shifted daily demand patterns due to the COVID-19 pandemic [2]. Additionally, it has become one of the most profitable options and is an energy resource that has recently decreased in cost. As a result, solar electricity generation has grown in residential, commercial, and utility-scale projects [3]. The future of PV generation will focus on optimizing hybrid systems [4] and improving the performance of each element of the system, as well as reducing their cost due to large-scale production [5]. Furthermore, the sector trend has been asking for low prices, and the competitive market has encouraged investment in solar PV technologies across the entire value chain, particularly in solar cells and modules, to improve efficiencies and reduce the levelized cost of energy (LCOE) [1]. As a result, PV power plants could grow almost sixfold over the next ten years, reaching a

**Citation:** Pérez-González, A.é.; Benítez-Montoya, N.; Jaramillo-Duque, Á.; Cano-Quintero, J.B. coverage path planning with Semantic Segmentation for UAV in PV Plants. *Appl. Sci.* **2021**, *11*, 12093. https:// doi.org/10.3390/app112412093

Academic Editors: Luis Hernández-Callejo, Maria del Carmen Alonso García and Sara Gallardo Saavedra

Received: 16 November 2021 Accepted: 15 December 2021 Published: 19 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/).

cumulative capacity of 2840 GW globally by 2030 and rising to 8519 GW by 2050, according to [6].

Thus, it is established that the number of PV plants and generated power have increased around the world. This implies specific technical challenges in their maintenance and operation (O and M) [7]. Some of these threats affect their production, which increases cost, and decreases profitability. The most common threats are the failures in the inverters and PV modules [8]. Through dirty equipment, the state of the environment, or manufacturing defects in the PV modules, PV plant energy generation can be curtailed by 31% in the worst cases [9–11].

One must consider that PV plants are commonly installed on roofs, rooftops, canopies, or facades for urban environments. Likewise, solar farms utilize rural environments, such as deserts, plains, and hills [7,12,13]. Depending on the location of the PV plant, the manual inspection tasks could be exhausting and take up to 8 h/MW, if the number of modules is considerable. The amount of inspection time increases for solar PV plants on rooftops or canopies, by virtue of aspects of their installation [14]. In addition, inspections to detect threats in the panels must be conducted by trained personnel. In some cases, problems occur in elevated installations, for which special training and certification is required. These jobs could put people and facilities at risk [15].

In recent decades, unmanned aerial vehicles (UAVs) have been increasingly used in inspection and patrol tasks [16–18]. UAV-based applications for PV plant inspections have many advantages in comparison with the manual inspection methods. The main advantages are flexibility, lower cost, larger area coverage, faster detection, higher precision, and the capacity to perform a superior and automatic inspections [11,17–21].

There are many approaches to performing an inspection with UAVs in PV modules. One of them used UAVs with a thermal imaging camera to take photos in the infrared spectrum to evaluate the UAVs parameters, such as height, speed, viewing angle, sun reflection, irradiance and temperature, all of which are necessary to perform defect inspection [22]. In a second approach, UAVs were implemented to inspect different solar PV plants, wherein analysis of the correlation between altitude and the pixel resolution was used to detect PV panel defects and features like shape, size, location, and color, among others, of a particular defect were also detected [7,23–26]. In a third approach, the authors proposed PV plant fault inspection with UAVs using image mosaics combined with orthophotography techniques to create a digital map; image mosaics were combined with multiple visible and/or infrared range images into a single mosaic image covering a large area [14,27,28]. Whereas the orthophotography technique is a vertical photograph that shows images of objects in true planimetric position [29]. Thus, these two techniques were integrated with previous works to achieve an advanced tool that allows monitoring and taking actions in the operation and maintenance (O and M) of PV plants [7]. In short, this tool has been used to detect defects and dust or dirt in PV modules [28]. Apart from this, some approaches have been developed for the detection of defects in PV modules using artificial vision techniques, machine learning, deep learning, and the integration therewith of the previous approaches [7,24,30–33]. Additionally, these approaches to planning the UAV's flight path were configured from a ground control station (GCS) program. Irrespective of these approaches, an important aspect is that the UAV should automatically follow the path to cover important points in PV plants. Many research efforts have been made to calculate the paths of and solve the waypoint planning problem for UAV inspection of solar PV plant applications [34–37], but none of them propose the coverage path planning (CPP) as a method to complete this task.

The CPP, given a region of interest (RoI) in a 2D environment, consists of calculating the path that passes through each one of the points that make up the desired environment and must be found considering the limitations of movement [38,39]. CPP is classified as a classical NP-hard problem in the field of computational complexity. These problems were initially analyzed for indoor environments with mobile robots. However, with the development of GPS, CPPs began to be used for missions with UAVs. Due to the environment in which the task is performed and the obstacles present, precise localization in the environment is an arduous task, which makes the CPP a difficult problem [40]. Additionally, it is classified as a motion planning subtopic in robotics, and has two approaches, heuristic and complete. In heuristic approaches, the robots follow a set of rules defining their behavior, but do not present a guarantee for successful full coverage. These guarantee using the cellular decomposition of the area, which involves space discretization into cells to simplify the coverage in each sub-region, unlike complete methods, which cannot afford such processing. Another important issue mentioned by the authors of reference [39] is the flight time to fully cover the area, which can be reduced using multiple robots and by reducing the number of turning maneuvers. Finally, the available RoI information is important; several approaches accept previous knowledge of the robot's respecting the search RoI (offline), while sensor-based approaches collect such information in real-time along the coverage (online) [41].

In the literature, CPP approaches are needed in several application areas, such as floor cleaning [42], agriculture [43,44], wildfire tracking [45], bush trimming [46], power line inspection [47], photogrammetry [39], visual inspection [48], and many more. Additionally, many surveys regarding CPP present several approaches and techniques for performing missions with, mostly, land vehicles [41,49,50]. The research interest in aerial robots (indoor and outdoor) has surely motivated the research of CPP [51]. This can be implemented in many UAVs platforms, such as fixed-wing, rotary-wing, and hybrid UAV (VTOL) [39]. Rotary-wing UAVs are inexpensive and have good maneuverability, and their small payload capacity limits the weight of on-board sensors and flight time. Hence, they are more suitable for CPP missions on a small scale. Additionally, the increasing usage of UAVs in applications with complicated missions has led to CPP methods being a very active research area for single and multiple UAVs, especially recently [39,52–54], As evidenced in a previous work [43], the classic taxonomy of coverage paths in UAVs are classified into no decomposition, exact cellular decomposition, and approximate cellular decomposition. The first performs the coverage with a single UAV, for which no decomposition technique is required, because the shape of the RoI has a non-complex geometry. The second divides the free space into simple, non-overlapping regions called cells. The union of all cells fills the free space completely. The cells in which there are no obstacles are easy to cover and can be covered by the robot with simple movements. The third is based on grids. They use a representation of the environment decomposed into a collection of uniformly squared cells [55], considering rectangular, concave, and convex polygons for RoIs. In addition to this, CPP performances were assessed with applied metrics according to [39]. The elemental approach most used to solve offline CPP problem sis the area decomposition into non-overlapping sub-regions [56], to determine the appropriate visiting sequence of each sub-region and to cover each decomposed region in a back-and-forth movement to secure a complete coverage path. As a result, the methods for obtaining complete coverage of an RoI are the exact and approximate cellular decomposition methods [57–59].

On the other hand, image processing helps to obtain a map of the robot or RoI. Robots, such as UAVs, need to know the RoI before commencing CPP [60], which represents where the PV plants are and can be determined in a process called boundary extraction [61,62]. Then, the deep learning (DL) image segmentation technique, also known as semantic segmentation [63], is achieved by applying deep convolutional neural networks (CNN), such as the U-Net network model [64,65] or the FCN model [66], which dramatically enhance the segmentation results. Once the segmentation is done and the mask is obtained or the RoI is identified, the GCS calculates the CPP that guides the UAV in the automatic plant inspection, during which it captures images of PV plants [67]. Most of the failures occur at the centimeter or millimeter level, and that poses a challenge for the inexpensive sensors available today [68]. Tests of coverage paths it can be conducted with drones in a real or simulated environment. The most viable option for this stage of the work is simulation, as verified by other research [69–72]. Simulation has been recognized as an important research tool; initially, simulation was an academic research tool, but with the

107

advancement of computers, simulation has reached new levels. It is a remarkable tool that guarantees support in design, planning, analysis, and making decisions in different areas of research and development. Of course, robotics, as modern technology, has not been an exception [73,74].

This work is focused on implementing the best strategy of coverage path planning (CPP) over PV plants with UAVs using semantic segmentation in a deep learning server to obtain the RoI. The experimental results were obtained by simulating the CPP methods and using UAVs, such as the 3DR Iris and Typhoon aerial robots.

The key contributions of this work are as follows:


This paper is structured as follows. In Section 2, necessary definitions and the techniques used to obtain the results are described. In Section 3, the three CPP methods implemented are compared to show relationships among the CPP width, number of maneuvers, and energy consumption, with the aim of finding the best CPP method to implement. Finally, in Section 4, some conclusions are given.

#### **2. Materials and Methods**

In this work, three PV power plants were selected that met the image requirements of no light distortions, non-complex geometry, and grouped panels. These plants are in different parts of the world. The first PV plant has an area of 35, 975 square meters, located in Brazil (−22.119948323621525, −51.44666423682603), known as Usina Solar Unioeste 1 [75]. The second PV plant has an area of 25, 056 square meters, located in Iran (34.0504329771808, 49.796635265177294), known as Arak power plant [37]. The third plant is in the United States (38.55989816199527, −121.42374978032115); the plant has 1344 square meters of area, and it is part of a PV plant located on the roof of the California State University Sacramento Library (CSUSL) [76].

These plants were subjected to a series of processing stages, as shown in Figure 1. The RoIs were obtained from Google Maps satellite images with a predefined altitude, and limits, from which the image (input image) of the desired PV system was obtained as the first stage, shown in Figure 1. In the second stage, the image was entered into a DL server that was developed in this work, taking into account the previous work [65]. The DL server was launched with TensorFlow [77] and Flask [78]. This stage obtains a mask (image segmentation) of the PV plant, as shown in Figure 1. In the third stage it was necessary to apply a series of OpenCV functions (post-processing), as referenced in Algorithm A1 in the Appendix A, to adjust the mask (output mask) to the PV plant area, as seen in Figure 1. Then, the output mask or RoI was introduced in the GUI interface, where the CPP method was selected to internally execute. Later, the path GPS positions (CPP computed) were sent to the UAV through MAVlink commands [79,80]. The UAV executed these commands in the Gazebo platform (CPP simulation) and, at the same time, the simulation data was fed back the GUI interface. Simultaneously all GPS points reached by the UAV were drawn in the GCS platform (CPP). In this stage, the trajectory was validated. Each stage is described in greater detail in the following sections.

**Figure 1.** Stages for CPP with semantic segmentation for UAV in PV plants.

#### *2.1. Deep Learning (DL) Server for Segmentation*

A deep learning (DL) server for segmentation was necessary to extract the RoIs from the Google Maps images. Additionally, the server was used to achieve this task automatically with the process called semantic segmentation, wherein each pixel is labeled with the class of its enclosing object [65,66]. In previous work, a convolutional neural network was proposed, wherein a public database was used; this data was prepared, resized for training, and assessed with two network structures. The U-net network had the best performance, in terms of metrics, in the semantic segmentation task [65]. Then a DL server was employed to perform image segmentation and obtain the RoI [81].

#### *2.2. Post Processing*

In this step, a set of OpenCV functions were applied with Python 3.7. For example, in the first function, morphological operators like "Erode" and "Dilate" were applied to the images, then the "FindCountours" function was applied to help extract their contours. The contour can be defined as a curve that joins all the continuous points at the boundary of the PV installation. So, the "ContourArea" function was then used to find the area of the previous contour. Following this pattern, the area was compared with 400 others to filter the bigger area and eliminate the little areas belonging to false positives. Then the "ApproxPolyDP" function was used to approximate a shape of the contour to another shape with fewer vertices. Subsequently, the "DrawContour" function was used to draw the resulting contour [82,83]. Finally, the "Erode" morphological operators were used again, to expand the known area and compensate for the limitations of the mask with regard to the CPP method and some faulty occurrences caused by the false positives of the DL server. The pseudocode of the openCV functions used is shown in Algorithm A1 in the Appendix A.

#### *2.3. 2D Coverage Path Planning Method in the GUI Interface*

In previous studies of CPPs, there were many existing methods from which to select to solve the CPP problem. In this work, three methods based on CPP were selected, considering the following criteria: time of execution, ease of implementation, and more, which were used to cover the RoI. The methods were selected according to [38,39].

The boustrophedon exact cellular decomposition (BECD), which was proposed by [84], was the first selected. The CPP method is delineated in Figure 2.

The second was grid-based spanning tree coverage (GBSTC), which works with cellular decomposition, first proposed by [85], and depicted in Figure 3.

**Figure 2.** The RoI and the path generated from the initial point (I) to the final point (F) by BECD.

The third method selected for this project was grid-based wavefront coverage (GBWC), first proposed by [86]. The method is illustrated in Figure 4. Each of these CPP methods is explained in more detail in the following sections.

(a). The boustrophedon exact cellular decomposition (BECD) Method: This method takes the robot's free space and obstacles and splits them into cells. These cells are covered by the robot using a back-and-forth pattern from the initial point to the final point, using maneuvers of 90 degrees to change direction from south to north or vice versa, as shown in Figure 2. This method improves the trapezoid decomposition technique, as it exploits the structure of the polygon to determine the start and end of an obstacle, and thus is able to divide the free space into a few cells that do not require a redundant step, and it permits the coverage of curved areas [84].

(b). Grid-Based Spanning Tree Coverage (GBSTC): This method is based on approximate cell decomposition and differs from the previous method in that the following postulates were considered. First, the method divides the space into grids of side *L*. Second, the robot only moves in perpendicular directions to the sides of the grid. Third, every grid is subdivided into four grids of side *L*/2. Finally, GBSTC discards space that is partly occupied by obstacles. Consequently, considering these previous postulates, the method consists of several stages: in the first stage, a graphic structure is defined, *S*(*N*, *E*), where *N* is nodes, defined as the central point of each grid, and *E* is edges, defined as the line segments connecting the centers of adjoining grids, as shown in Figure 3a. In the second stage, the method builds a spanning tree for S, and employs this tree to plan a cover path as follows. Starting in grid *I* with a sub-grid of side *L*/2, the robot begins by travelling between adjoining sub-grids along a path that moves around the spanning tree at a constant distance and in a counterclockwise direction, finishing when the initial sub-grid, *I*, is found again, which means it is also in the final point, *F* [85]. An example of this method is illustrated in Figure 3b. The approximation depends on the side length, *L*, of the grid.

**Figure 3.** (**a**) Approximate cell decomposition in grids, and sub grids. (**b**), coverage path generated with the GBSTC method.

(c). grid-based wavefront coverage (GBWC): The first grid-based method proposed for CPP, GBWC is an offline method that uses a grid representation, in addition to applying a full CPP method. The method requires an initial grid, *I*, and a final grid, *F*. A distance transformation that propagates a wavefront from the final to the initial point is used to assign a specific number to each item on the grid. That is, the method first assigns a zero to the final item, and then a one to all its surrounding grids. Then all unmarked grids adjoining those marked one are numbered two. The process is incrementally repeated until the wavefront reaches the initial grid [86], as illustrated in Figure 4a.

**Figure 4.** (**a**) Wavefront distance transforms for the selection of the initial position (I) and final position (F). (**b**) Coverage path generated using the wavefront distance transforms with the selection of the initial position (I).

Once the distance transformation is determined, a coverage path can be found by starting at the initial grid, *I*, and selecting the adjoining grid with the highest number that has not been explored. If two or more are unexplored, and adjoining grids share the same number, one of them is randomly selected, as shown in Figure 4b.

#### 2.3.1. Metrics

The metrics evaluate the performance of the three CPP methods. Such assessment can be performed by considering five commonly used metrics to evaluate the effectiveness of the proposed CPP methods both theoretically and in simulation (dynamically) [7,41,43,62]. These five metrics are covered path length, flight time, energy consumption, redundancy of points traveled, and percentage of coverage of the total area. Each of these metrics is described below.

The *Lct* metric (covered path length) is the length of the entire path covered by the UAV from the initial to the final point. For a trajectory in a 2D plane composed of *n* points, assuming the initial point as (*x*1, *f*(*x*1)) and the end point as (*xn*, *f*(*xn*)), the *Lct* can be computed as shown in Equation (1):

$$L\_{ct} = \sum\_{i=1}^{n-1} \sqrt{(\mathbf{x}\_i + \mathbf{1} - \mathbf{x}\_i)^2 + (f(\mathbf{x}\_1 + \mathbf{1}) - f(\mathbf{x}\_i))^2} \tag{1}$$

where (*x<sup>i</sup>* , *f*(*xi*)), in which *i* = 1, 2, . . . , *n*, representing all the *n* points of the UAV flight path in 2D coordinates. More details can be found in [87].

The flight time metric of the PX4 SITL Gazebo model is the time required to travel the total flight path (takeoff, path travel, and landing) with a dynamic speed that considers the inertia and the variation in speed due to UAV turns angles. These data are collected through sensors in the Gazebo plugins [88].

The redundancy of points traveled, *R*%, corresponds to the number of points that are visited more than once from the total number of points that the path contains, Equation (2).

$$R\_{\%} = \frac{P\_{pc}}{P\_{vmo}} \tag{2}$$

where *Ppc*, and *Pvmo* correspond to the number of points the trajectory contains, and the number of points visited more than once.

The percentage of coverage of the total area, *C*%, measures the number of effectively covered points within the total number thereof by the points the area to be covered contains, given by Equation (3).

$$\mathcal{C}\_{\%} = \frac{P\_v}{P\_T} \tag{3}$$

where *Pv*, and *P<sup>T</sup>* correspond, respectively, to the total number of points visited and the total number of points in the area.

Moreover, the number of maneuvers metric, which is the number of turning maneuvers the UAV performs on a path, is often used as the main performance metric in coverage [89,90].

The energy consumption metric is computed from the voltage and current data from the power module of the PX4 SITL model [91]. Its value depends on parameters, such as the CPP width between lines and the speed and the height of the UAV at the time of implementing the CPP, which were configured in the interface and were simulated in Gazebo.

To validate the results of the methods described above, the BECD, GBSTC, and GBWC methods were implemented in two UAVs, simulated in Gazebo, and their performances were assessed by the metrics presented above [39]. The next section describes the results and compares the models in detail.

#### *2.4. Simulation and Validation Platform*

Based on [92], the Gazebo platform was selected to execute the simulation experiments, as it has extensive documentation on its webpage. In addition, it is the most mentioned and used simulation platform in previous work [73,74] that implemented path planning or CPP, used UAV sensors, or deployed several UAVs [92]. Gazebo also allows the modeling of commercial UAVs using the PX4 autopilot software [91], as shown in Figure 5. In addition, this figure shows one of the experiments conducted with the Typhoon UAV, flown over the CSUSL plant. The UAV sonar sensor, represented by blue lines, is also shown, as is the UAV camera, in the box at the upper right.

**Figure 5.** UAV over A PV plant, simulated in Gazebo.

The integration of Gazebo [92], PX4 [91], Python, and the CPP methods was implemented using ROS (Robot Operating System) as middleware [93]. This tool allows communication among nodes. The nodes are processes, and each node has a task associated with it, such as sending a Mavlink command to control the UAV trajectory and permitting reading messages from the UAV to discern its status in flight, using a simulation mode referred to as software in the loop (SITL). This simulator provides the ability to run different vehicles, such as a plane, copter, or rover, without a need for any microcontrollers or hardware [94]. In addition, two PX4 autopilot rotary wing UAVs, the 3DR Iris and Typhoon UAV, were chosen because they have good maneuverability and are more suitable for small-scale CPP missions. Furthermore, these UAVs have been widely used in other research [67,70,95].

The GCS software (QGroundControl) was selected to validate the CPP calculated for each UAV and PV plant obtained, because it is the most compatible tool with the PX4 autopilot. It is also recommended on the PX4 webpage. On other hand, another compatible tool, namely a GUI, was designed using Qt to modify and vary the parameters and to convert the way points from RoI pixels to geo-referenced points [96].

#### **3. Results and Discussion**

The procedure described above—DL-server, post-processing, CPP, and metrics evaluation—are represented in Figure 1, was applied to three different PV plant images. The results obtained are described in the following sections.

#### *3.1. Results with Deep Server, and OpenCV Functions*

OpenCV functions and a DL server were combined to accurately extract the mask and to trade off the DL server results of the PV plant area. Then, the CPP methods used the mask as a region of interest (RoI).

The stages to obtain the results of the RoI in the three images are shown in Figures 6–8. In the first stage, a high-resolution image from Google Maps, of predefined height and width, was obtained, and used as an input image in the DL server, depicted in Figure 6a. In the second stage, the output image of the DL server was the mask, shown in Figure 6b. In the third stage, the opening function was used in the mask, Figure 6c. In the fourth stage, the RoI was obtained using the draw contour method, as seen in Figure 6d. Finally, the RoI was blended with the input image with the aim of comparing the results in Figure 6e. The results were satisfactory and can be adapted, depending on the environment.

**Figure 6.** Steps of boundary extraction by the DL server and OpenCV functions in the Unioeste 1 PV plant.

**Figure 7.** Steps of boundary extraction by DL server, and OpenCV functions in Arak PV plant.

**Figure 8.** Steps of boundary extraction by DL server, and OpenCV functions in CSUSL PV plant.

#### *3.2. Results of the CPP Method*

The CCP method results were obtained from six experiments that implemented the five metrics previously described. Each experiment was conducted by selecting a PV plant and a UAV, then choosing the CPP method and the width, speed, and height of the UAV over the flight path. All of this was done in the implemented GUI interface. The six tests were conducted by varying the CPP width parameter, while other parameters were not varied, and the same experimental conditions were maintained throughout.

For these experiments two UAVs were selected, the Typhoon [97] and 3DR Iris [98], as mentioned in Section 2.4. Three simulated PV plants were also chosen, as highlighted in Section 2. The metrics referenced in Section 2.3.1 were assessed for each test. Battery consumption was obtained from the GCS software with the SITL parameter activated. The rest of the data was collected from the Gazebo simulations. The experiments are explained in detail in the following sections, and some results are discussed.

#### 3.2.1. The Three First Experiments with a Typhoon UAV

The three first CPP experiments were conducted with a Typhoon UAV, varying the CPP width between 0 and 20 m, depending on the PV plant being covered. For example, in the first and second PV plants, Unioeste 1, and Arak, respectively, the CPP widths were varied between 5 and 20 m. In the third CSUSL PV plant, the width was varied between 1 and 8 m, due to the method's restrictions on running in small areas with a large width. Then, in each of these experiments, the resulting metrics were annotated in Tables, as shown in Appendix B.

The logged metrics were utilized to draw a clustered columns and lines chart for experiments 1, 2 and 3, performed with the Typhoon UAV, to highlight the most important information and to determine their correlations with each other. In these graphs the vertical axis was scaled logarithmically and the horizontal axis was labeled with the type of CPP method in use; in this way, it is possible to visualize the correlation of flight time, number of maneuvers, and path length covered with the energy consumption by each CPP method, as shown in Figures 9–11. In addition, the relationship between the redundancy metric and the CPP width is identified, showing that this metric is higher in BECD, but does not affect energy consumption, which is the primary metric of interest.

**Figure 9.** Performance metrics of Typhoon UAV on Unioeste PV plant (Experiment 1).

**Figure 10.** Performance metrics of Typhoon UAV on Arak PV plant (Experiment 2).

**Figure 11.** Performance metrics of Typhoon UAV on CSUSL PV plant (Experiment 3).

The UAVs' energy consumption, in conjunction with the CPP widths, were obtained from the energy consumption results logs, then a new Table 1 was composed of three large columns, corresponding to all the assessed PV plants, and each PV plant column is made up of three CPPs with their respective values of consumed energy. The rows of the table contain the CPP widths in increasing order from top to bottom, as shown in Table 1. Where it can be observed, some parameters, such as the energy consumption with regard to CPP width, were due to the size of the RoI. Additionally, it can be observed that the greater the CPP width, the lower the energy consumption, as seen in the columns of Table 1. For example, for the BECD CPP in Unioeste 1, with a width of 5 m, the percentage of energy consumed was 98%, while, for a width of 20, the energy consumed was 29%. It is also observed that the larger the PV plant, the UAV consumed all its energy in CPPs with narrow widths; on the other hand, if the plant is small, with the same width, the UAV does not have energy consumption problems. As can be seen in the row with a width of 8, where, for the Unioeste 1 and Arak PV plants, a lot of energy was consumed, between 88% and 52%, unlike the CSUSL plant, in which where the energy consumed was very little, between 7% and 5%. In the experiments for which energy consumption, with respect to CPP width, cannot be observed, (*N*/*A*) was used to annotate these results, which happened when the CPP width was very large with respect to the RoI, as this does not allow the generation of the route, or when the CPP width was very small with regard to the RoI causing a flight path in which the UAV consumes all its energy. In short, the UAV used can be undersized or oversized with respect to its intended PV plant.

The graphs in Figure 12 were constructed from Table 1 by quintic polynomial interpolations, which requires six data points to form a curve that passes through all given data points [99], where the abscissa axis is the CPP width and the ordinate axis is energy consumption. For the experiment conducted at the PV plant Unioeste 1, the graph in Figure 12a shows that the BECD method had the lowest energy consumption when the CPP width was in the range of 5 to 16 m, while GBSTC and GBWC had the lower energy consumption when the CPP width was in the range of 16 to 20 m. For the experiment tested at the Arak PV plant, the graph in Figure 12b shows that the BECD method had the lowest energy consumption when the CPP width was in the range of 5 to 10 m, while the GBSTC and GBWC had similar energy consumption when the CPP width was in the range of 10 to 15 m, and GBWC had the lowest energy consumption when the CPP width was in the range of 15 and 20 m. For tests conducted at the CSUSL PV plant, the graph in Figure 12c shows that the BECD method had the lowest energy consumption when the CPP width was in the range of 1 to 5.5 m, while the method GBSTC had the lowest energy consumption when the CPP width was in the range of 5.5 to 8 m. All plants were recreated in a simulation environment with their real dimensions.


**Table 1.** Comparison of three CPP methods with regard to CPP width, for a simulated Typhoon UAV at three PV plants.

**Figure 12.** Performance test of Typhoon UAV varying the CPP width.

#### 3.2.2. The Last Three Experiments with 3DR Iris UAV

The last three CPP experiments were conducted with a 3DR Iris UAV, varying the CPP width between 1 and 20 m, depending on the PV plant covered. For example, in the first and second PV plants, Unioeste 1 and Arak, respectively, the CPP widths were varied between 8 and 20 m; in the third PV plant, CSUSL, the width was varied between 1 and 8 m, due to the restrictions of the method requiring it be run in small areas with a large width. Then, in each of these experiments, the resulting metrics were annotated in tables, as shown in Appendix C.

As in the previous chapter, the recorded metrics were used to create bar and line graphs for Experiments 4, 5 and 6, which are those conducted with the 3DR Iris UAV, to highlight the most important information and to determine correlations among the metrics. In these graphs, the vertical and horizontal axes were scaled logarithmically and labeled in the same manner as above, with the aim of visualizing the correlations between flight time, number of maneuvers and path length covered, with regard to energy consumption by each CPP method, as shown in Figure 13–15 . Furthermore, the relationship between the redundancy metric and CPP width is identified, showing behavior similar to the previous UAV's data.

**Figure 13.** Performance metrics of 3DR Iris UAV on Unioeste PV plants (Experiment 4).

**Figure 14.** Performance metrics of 3DR Iris UAV on Arak PV plants (Experiment 5).

**Figure 15.** Performance metrics of 3DR Iris UAV on CSUSL PV plants (Experiment 6).

The UAV energy consumption and CPP widths were obtained from the energy consumption results logs, then a new Table 2 was established, in which three large columns correspond to the PV plants, where each PV plant column contains three CPPs with their respective values of consumed energy, and the rows of the table contain the CPP widths in an increasing direction from top to bottom. Here it can be observed, as in the last table, that

some parameters' values, such as the energy consumption by CPP width, are due to the size of the RoI. It can also be observed that the greater the CPP width, the lower the energy consumption, as seen in the columns of Table 2. For example, using BECD at Unioeste 1 with a width of 10 m, the percentage of energy consumed was 100%, while for a width of 20 m, the energy consumed was 48%. It is also observed that, for larger PV plants, the UAV consumed all its energy in a CPP with a small width, whereas, for small plants of the same width, the UAV did not have energy waste problems. As can be seen in the row corresponding to a CPP width of 8 m, at the Unioeste 1 and Arak plants a lot of energy was consumed, between 100%, and 88%, unlike at the CSUSL plant, where the energy consumed was very little, between 9%, and 8%. Similarly, in the previous experiments, wherein the UAV's energy was exhausted and the CPP width did not allow the generation of the route or when the CPP width was large with regard to the RoI, the results were scored with (*N*/*A*). To summarize, the UAV used is suitable for the last PV plant.

On the other hand, the 3DR Iris UAV, with its design, size, autonomy, and performance is suitable when the PV plant to be inspected is small, of an approximate size of 5000 square meters or less, such as the CSUSL PV plant, at 1344 square meters. Owing to such PV plant sizes, any coverage method used in this work can be used to cover an area with a CPP width of 1 meter, while minimizing the metrics of path length covered, flight time and maneuverability, to obtain lower energy consumption, as shown in the Table 2. In addition, this UAV did not perform well at large or medium-sized plants, since, to cover them using the CPP method, the width that must be provided is greater than 8 m; therefore, though the UAV's inspection of the PV plant had no high-resolution cameras on board, they were needed; without them, the inspection cannot be guaranteed.


**Table 2.** Comparison of three CPP methods with regard to CPP width, simulated in a 3DR Iris UAV on three PV plants.

The graphs in Figure 16 were constructed from Table 2 using quintic polynomial interpolations, wherein the abscissa axis is CPP width and the ordinate axis is energy consumption.

Additionally, the graph in Figure 16a shows that the BECD method had the lowest energy consumption; when the CPP width spanned the entire test range, the GBSTC and GBWC methods showed a higher energy consumption at the Unioeste 1 PV plant. The graph in Figure 16b shows that the BECD method had the lowest energy consumption when the CPP width was in the range of 8 to 10 m, whereas BECD and GBWC had similar energy consumption, leaving the GBWC method with the lowest consumption, when the

CPP width was in the range of 10 to 15 m. On the other hand, the GBWC method had the lowest energy consumption when the CPP width was in the range of 15 and 20 m; this test was conducted at the Arak PV plant. Finally, the graph in Figure 16c shows that the BECD method had the lowest energy consumption when the CPP width was in the range of 1 to 3 m, and the GBWC method had the lowest consumed energy when the CPP width was in the range of 3 to 8 m. The test was performed at CSUSL's PV plant.

**Figure 16.** Performance test of UAV 3DR Iris varying the CPP width.

Summary Tables of the results of the metrics for each of the tests with the Typhoon UAV and 3DR Iris UAV are shown. All files, and logs for the experiments are available on GitHub at [100].

#### *3.3. Discussion*

The proposed strategy allows a semi-automatic and faster solution for achieving effective results when inspecting PV plants in geometrically simple areas, with certain limitations, although some results obtained in this work are theoretical results, such as from CPP width, for example; in practice, a camera of 14 megapixels is not adequate to inspect a PV plant with a CPP width of 20 m [37,101].

The results obtained in this work indicate that the most adequate method is BECD for a specific range of CPP widths, although it also shows adequate performance for various CPP widths in some RoIs. The GBWC showed a good performance when using a CPP width greater than 7 m in some RoIs, and all showed a relationship with the area to be covered, as shown in Tables 1 and 2.

An analysis of the metrics used in this work, such as redundancy, *R*%, which does not represent a significant factor when comparing these three CPP methods, showed that, although the BECD was the method with the highest redundancy, it was also the method with the lowest consumption of energy. On the other hand, the other metrics, such as the *Lct* metric, flight time metric, and the number of maneuvers, are directly related to energy consumption, as can be seen in Appendix B and C. It also helps to conclude that the BECD method is more suitable for widths in a range between 0 and 7 m, due to the lower number of maneuvers in these ranges, as other authors have also mentioned [102,103]. The percentage of coverage of the total area, *C*%, always showed that coverage was total

The *Lct* metric, flight time and number of maneuvers were directly related to the consumption of energy for all the experiments performed, as seen in Figures 9–11 and 13–15, and as confirmed by other authors [89,90], in addition to helping to reinforce that the most appropriate CPP for certain ranges is the BECD.

The BECD method is the best method among the three methods tested, in a specific width range, from 0 to 7 m, for all the RoIs tested, which means that for a 10 Megapixel camera with a horizontal field of view of 7 m, the CCP method's good images, in terms of what to inspect of a PV plant can be obtained [23].

The implemented BECD method divided the RoI into small regions, and then, over these regions, it implemented the round-trip coverage pattern, called boustrophedon [38,39]. This pattern allows spending less energy consumption with small widths, due to its low number of maneuvers compared with the other methods, as shown in Appendix B and C.

On the other hand, the GBWC method allowed lower energy consumption, according to the results of the simulations, in widths greater than a certain value, but depending on RoI size; this method also allows an approximate coverage outside the RoI, a characteristic that is very important in this type of application, and that, perhaps, is not very attractive for terrestrial robotics, from which this type of method originates [86].

A RoI with many concave points is a great challenge for performance metrics, since they increase with the distance of travel, and also increase the number of maneuvers and therefore increase energy consumption, according to these results in the following works [67,103], a more detailed analysis on these characteristics was made.

In future works, the methods (BECD, GBWC) could be implemented in UAVs with characteristics similar to those used in this work; these characteristics are shown in Tables A1 and A5. Additionally, with these UAVs, an inspection could be conducted in at least one PV plant. Where one can think about the implementation of an expert system that selects the coverage path between the two types of methods (BECD, GBWC), according to the CPP width size, required for a camera with a certain resolution, and focal length.

Finally, the results obtained with the BECD, and GBWC methods differ from the results obtained by other authors [34–37]. Who did not implement CPPs to solve these types of problems; they used other types of solutions that have restrictions when inspecting PV plants with UAVs. On the contrary, this work considers the CPPs, and obtained interesting results for future real implementation. The proposed CPP will increase the possibility of using inexpensive UAV systems for the inspection of PV systems on roofs of houses, and commercial buildings, and also, of the use of CPP with small widths to complete inspections at centimeter scales of the panel where the flaws can be better seen.

#### **4. Conclusions**

In this work, a method for implementing CPP in UAV for PV plant inspection was presented. The method consisted of a series of steps, one of these was the deployment of a DL-based U-net model to establish a DL service system from which to extract the limits of PV plants by extracting the boundaries of PV plants from an image. To summarize, the method was accurate, and fast without depending on the image, with low request latency and response.

This experiment focused on three novel path planning methods in the PV inspection missions in order to find the best path for covering each of the three PV power plants with less energy consumption. A GUI interface was used to order the UAV's maneuvers in the inspection of the simulated PV plants. The results of each CPP method in the simulation were compared. The best CPPs was the BECD, for a range of CPP widths of 0 to 7 m. These path planning algorithms can be performed by any multirotor UAV that receives Mavlink commands, can carry a camera sensor, and transmit real-time video to GCS.

Performance on the CPP tasks was measured using two different types of flying robots, a Typhoon UAV and a 3DR Iris UAV. It was demonstrated that the Typhoon UAV (or one with similar characteristics) is better suited for large or medium-sized PV plants (such as those in deserts, plains, and hills); instead, the drone-like 3DR Iris UAV is more suitable for small-sized PV generation (such as on roofs and rooftops, canopies, and facades). The proposed strategy allows such comparisons to be made and enables the selection of the most suitable UAV for each type of installation.

The values obtained for the metrics collected from each of the tests show a correlation between covered path length, flight time, number of maneuvers as regards energy consumption, and ensuring that the CPPs implemented in UAVs to inspect photovoltaic plants could be similar when implemented in real plants. The results also help to predict the energy consumption of a given UAV when performing a plant inspection.

**Author Contributions:** Conceptualization, A.é.P.-G., N.B.-M., Á.J.-D. and J.B.C.-Q.; methodology, A.é.P.-G.; software, A.é.P.-G., N.B.-M.; validation, A.é.P.-G.; formal analysis, A.é.P.-G.; investigation, A.é.P.-G., N.B.-M.; resources, A.é.P.-G., Á.J.-D. and J.B.C.-Q.; data curation, A.é.P.-G.; writing–original draft preparation, A.é.P.-G.; writing–review and editing, A.é.P.-G., N.B.-M., Á.J.-D. and J.B.C.-Q.; visualization, A.é.P.-G.; supervision, Á.J.-D. and J.B.C.-Q.; project administration, A.é.P.-G., Á.J.-D. and J.B.C.-Q.; funding acquisition, A.é.P.-G., Á.J.-D. and J.B.C.-Q. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Colombia Scientific Program within the framework of the so-called Ecosistema Científico (Contract No. FP44842-218-2018).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The models used in the computational experiment are available at GitHub in [100].

**Acknowledgments:** The authors gratefully acknowledge the support from the Colombia Scientific Program within the framework of the call Ecosistema Científico (Contract No. FP44842-218-2018). The authors also want to acknowledge Universidad de Antioquia for its support through the project "estrategia de sostenibilidad".

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A. Algorithm**

In this step, a series of OpenCV functions, in python 3.7 and in jupyter-lab, were applied.


#### **Appendix B. Tables of Typhoon UAV**

**Table A1.** Yuneec Typhoon UAV technical specification.


**Table A2.** Experiment 1 with three CPP methods implemented using Typhoon UAV over Unioeste 1 PV Plant.

**Table A3.** Experiment 2 with three CPP methods implemented using Typhoon UAV over Arak PV Plant.



**Table A4.** Experiment 3 with three CPP methods implemented using Typhoon UAV over CSUSL PV Plant.

#### **Appendix C. Tables of 3DR Iris UAV**

**Table A5.** 3D Robotics UAV 3DR Iris technical specification.



**Table A6.** Experiment 4 with three CPP methods implemented using UAV 3DR Iris over Unioeste 1 PV Plant.

**Table A7.** Experiment 5 with three CPP methods implemented using the UAV 3DR Iris over Arak PV Plant.



**Table A8.** Experiment 6 with three CPP methods implemented using the UAV 3DR Iris over CSUSL PV Plant.

#### **References**

	- 93. Quigley, M.; Conley, K.; Gerkey, B.; Faust, J.; Foote, T.; Leibs, J.; Wheeler, R.; Ng, A.Y. ROS: An open-source Robot Operating System. In Proceedings of the ICRA Workshop on Open Source Software, Kobe, Japan, 12–17 May 2009; p. 5.
	- 94. Silano, G.; Iannelli, L. CrazyS: A Software-in-the-Loop Simulation Platform for the Crazyflie 2.0 Nano-Quadcopter. In *Robot Operating System (ROS): The Complete Reference (Volume 4)*; Koubaa, A., Ed.; Springer International Publishing: Cham, Swotzerland, 2020; pp. 81–115. [CrossRef]
	- 95. Salamh, F.E.; Karabiyik, U.; Rogers, M.K. RPAS Forensic Validation Analysis Towards a Technical Investigation Process: A Case Study of Yuneec Typhoon H. *Sensors* **2019**, *19*, 3246. [CrossRef]
	- 96. QGroundControl Ground Control Station. 2021. Available online: https://github.com/mavlink/qgroundcontrol (accessed on 31 October 2021).
	- 97. Yuneec Typhoon H-Yuneec Futurhobby. Available online: https://yuneec-futurhobby.com/yuneec-typhoon-h-pro-realsense (accessed on 28 October 2021).
