Energy Management System for the Optimal Operation of PV Generators in Distribution Systems Using the Antlion Optimizer: A Colombian Urban and Rural Case Study

Abstract

This paper presents an Energy Management System (EMS) for solving the problem regarding the optimal daily operation of Photovoltaic (PV) distributed generators in Alternate Current (AC) distribution grids. To this effect, a nonlinear programming problem (NLP) was formulated which considered the improvement of economic (investment and maintenance costs), technical (energy losses), and environmental (${CO}_2$ emission) grid indices as objective functions, considering all technical and operating constraints for the operation of AC networks with the presence of PV sources. To solve this mathematical formulation, a master–slave methodology was implemented, whose master stage employed the antlion optimizer to find the power dispatch of PV sources in each period of time considered (24 h). In the slave stage, an hourly power flow based on the successive approximations method was used in order to obtain the values of the objective functions and constraints associated with each possible PV power configuration proposed by the master stage. To evaluate the effectiveness and robustness of the proposed methodology, two test scenarios were used, which included three installed PV sources in an urban and a rural network, considering the PV power generation and demand located reported for Medellín and Capurganá, respectively. These systems correspond to connected and standalone grids located in two different regions of Colombia. Furthermore, the proposed methodology was compared with three optimization methodologies reported in the literature: the Chu and Beasley genetic algorithm, the particle swarm optimization algorithm, and the vortex search optimization algorithm. Simulation results were obtained via the MATLAB software for both test scenarios with all the optimization methodologies. It was demonstrated that the proposed methodology yields the best results in terms of solution quality and repeatability, with shorter processing times.

1. Introduction

1.1. General Context

In recent years, the growing dependence of human beings on electrical systems has led to the use and exploitation of fossil fuels in order to meet the demand for electrical energy, with negative impacts on the environment, since polluting gases are emitted into the atmosphere [1,2,3]. In order to overcome this problem, various governments worldwide have promoted the integration and use of power generation from renewable sources. The most widely used renewable energy sources worldwide are solar panels and wind turbines, as they are considered to provide clean and unlimited energy (in the scale of human energy consumption), in addition to the fact that their acquisition costs have been reduced thanks to the different technological advances made in the last decades [4,5,6].

In the Colombian context, approximately 68.3% of electricity is generated from hydroelectric plants that depend directly on the country’s weather conditions [7]. It is for this reason that thermal power plants play an important role, as not only do they support hydroelectric plants but also cover 30.7% of the power demanded by users. However, these plants use fossil fuels such as natural gas (13.3%), coal (9.6%), and diesel (7.8%), which contribute to the emission of polluting gases [7]. Therefore, the integration of energy generation sources based on renewable sources for electrical systems has been promoted in the national territory (Colombia) by means of regulations and legislation, e.g., Law 1715 of 2014 and CREG Resolutions 030 of 2018 and 068 of 2020.

On the other hand, due to the geographical location of the country, it has been possible to use, develop, and generalize the use of renewable energies based on solar resources [8]. This has allowed Colombia to start developing large-scale projects related to photovoltaic (PV) solar energy [8]. Although the installed capacity is currently far from the maximum usable levels, this has made it possible to propose incentives and updates in the form of laws, standards, regulations, and policies for the rational and efficient use of renewable energies [8]. However, due to the variation and uncertainty of the solar resource, which is caused by weather conditions that depend on the location of the electrical system, as well as on the period of the year [6,9], the operation of PV generators in traditional electrical networks and in electrical distribution grids poses many technical and operations challenges, as an incorrect operation can entail stability and reliability issues in the electrical system [10]. To deal with this problem, it is necessary to implement Energy Management Systems (EMS). One of the main objectives of these systems is to determine the optimal power dispatch of energy sources based on renewable resources to economically feed system loads while ensuring a high-quality service that is reliable and environmentally sustainable [11,12,13]. An EMS can operate in two different ways: in centralized and decentralized operation schemes [11]. In this research, a decentralized EMS is implemented which seeks to provide PV generators with a greater degree of freedom, thus achieving a correct balance between generation and demand for all periods of operation of the distribution system.

1.2. Motivation

In the Colombian context, the isolated or non-interconnected zones (ZNI) correspond to 52% of the national territory, where, due to their difficult access, are energetically isolated from the rest of the national territory, and the costs of providing the electricity service are high when compared to urban areas that are connected to the National Interconnected System (SIN). This affects more than 2 million people (out of which more than 1.2 million do not have access to electricity), causing energy generation to be produced mainly through diesel generation plants [14]. It is important to note that the diesel fuel used in these remote areas is efficient and easy to obtain in the market, but its use produces great amounts of polluting gases that contribute to global warming. Thereupon, the main motivation of this research is to propose an EMS that allows correctly and intelligently managing the power injected by PV generators, which will allow distribution grids to deliver a high-quality electrical service that is as economically as possible and has a low environmental impact to each one of their users not only in urban areas but also in rural or isolated areas.

1.3. Literature Review

Energy management in a distribution grid can be oriented toward meeting technical, economic, or environmental criteria [15]. Thus, in the specialized literature, different techniques and approaches have been proposed for the implementation of EMS. Some of the most recent publications are presented below.

The authors of [16] propose a second-order conic programming model to optimally manage an interconnected distribution grid. The main objective of this document is to minimize CO${}_2$ emissions and the energy loss costs for one day of operation. The results obtained show the applicability and efficiency of the proposed methodology. In [17], an EMS for dispersed generators dispatch is proposed for distribution grids operating with direct current under a specified demand condition. The main objective is to reduce power losses in different penetration scenarios using dispersed generators. To this effect, a master–slave methodology is presented, which is based on the black hole algorithm in conjunction with the Gauss–Seidel method. The authors of this document validate the applicability of the proposed methodology by comparing the results obtained with three other metaheuristic algorithms and the GAMS software. However, a statistical analysis aimed at determining the repeatability of the solutions was not performed nor were the required processing times analyzed. Likewise, the minimization of power losses through dispersed generators dispatch in DC distribution grids under a given demand condition has also been solved by implementing sequential quadratic programming models [18] and convex quadratic programming models [19]. In these works, the authors demonstrate the efficiency of the proposed approaches by comparing their results with different methodologies reported in the literature that are mainly based on metaheuristic optimizers. The authors of [20] address the distribution grid management problem through a nonlinear programming model whose objective is to minimize CO${}_2$ emissions and the costs of energy losses for one day of operation. In this case, the distribution grid is made up of photovoltaic panels, wind turbines, and energy storage devices. This model is solved by using the specialized GAMS software, with which it is possible to observe significant reductions in polluting gases and energy losses.

To better illustrate the publications associated with EMS in distribution grids,

Table 1. Main reports involving EMS applied to electrical distribution grids.

Reference	Methodology	Indicator	Computation Time	Statistical Analysis	Comparison with Other Methodologies
[21]	HOMER pro	Technical–economic	No	No	Yes
[22]	MATLAB Optimization Toolbox	Technical	No	No	No
[17]	Master–Slave	Technical	No	No	Yes
[18]	Quadratic Programming	Technical	Yes	No	Yes
[19]	MATLAB CVX	Technical	No	No	Yes
[23]	Master–Slave	Technical	No	No	Yes
[24]	Fuzzy Logic	Environmental	No	Yes	Yes
[16]	MATLAB CVX	Environmental–economic	No	No	No
[20]	GAMS	Environmental–economic	No	No	No
[25]	Genetic Algorithm	Economic	No	No	No
[26]	MATLAB linprogParticle Swarm	Economic	No	No	No
[27]	Quadratic Programming	Economic	No	No	No
[28]	Universal Generator Function	Technical–economic	No	No	Yes
[29]	Online EMS	Technical–economic	No	No	No
[30]	Particle Swarm	Economic	No	No	No
[31]	Metaheuristic Algorithms	Economic	No	Yes	Yes
[32]	Artificial Fish Swarm	Economic	No	No	No
[33]	Particle Swarm	Technical–economic	No	No	No
[34]	Dynamic Programming	Economic	No	No	Yes
[35]	Master–Slave	Economic	Yes	Yes	Yes
[36]	Genetic Algorithm	Technical–economic	No	Yes	Yes
[37]	HOMER pro	Technical–economic	No	No	No
[38]	SAM	Technical–economic	No	No	Yes
[39]	SAM	Technical–economic	No	No	Yes

summarizes the main energy management strategies reported in the literature.

As seen in Table 1, the EMS seeks the correct operation of power generation devices that may be present in an MG, so that the owner and/or operator of the network can obtain greater benefits both in technical, economic, and environmental aspects while observing the technical-operative conditions of the system. It was also possible to identify a greater interest in methodologies that seek to improve the economic indicators of distribution grids, with the improvement of technical and environmental indicators being currently under development. Likewise, a high use of specialized software can be appreciated, as is the case of HOMER pro, SAM, GAMS, and CVX, in order to implement EMS. In the same way, statistical analysis and computation times are not widely adopted criteria nor is the use of comparison methodologies, which makes it difficult to identify the effectiveness and speed of the proposed solution strategies. In addition, the vast majority of the methodologies presented in this review of the state of the art do not consider the telescopic characteristic of radial networks, which does not allow representing what actually happens with existing electrical networks within a distributed generation environment.

1.4. Contribution and Scope

This research proposes the implementation of an EMS designed upon the basis of a master–slave methodology, which seeks to improve the technical, economic, and environmental indicators associated with the operation of electrical distribution grids. Additionally, the mathematical model includes the set of constraints that represent the operation of a distribution grid in a PV generation environment, i.e., the balance of active and reactive power, voltage regulation, the thermal limit of the conductors, the variations in the solar resource, and typical demand behaviors. Based on the review presented above, this study makes the following contributions:

i.

The use and presentation of a detailed mathematical formulation that contemplates the thermal limit of distribution grid conductors, which allows representing the operation of PV generators in telescopic radial networks with a higher degree of realism. The objective function of this mathematical model is the minimization of the technical, economic, and environmental indicators, and the set of constraints manages to capture the behavior of a distribution grid in a PV generation environment.

ii.

The implementation of a new master–slave methodology that allows designing an EMS for the optimal dispatch of PV generators. The master stage uses the antlion optimizer (ALO) to define the power injection of the PV generators in the MG. Regarding the slave stage, the successive approximations power flow method is used to evaluate the technical, economic, and environmental indicators associated with the operation of the system.

iii.

A methodology that finds the global optimal solution for a complex optimization problem from a dimensional perspective, such as the PV generator operation problem in distribution grids, thus achieving the best results in terms of solution quality and repeatability.

Due to the fact that Colombia is a country located between the tropics of Cancer and Capricorn, its energy consumption habits are not influenced by seasonal changes, so the resource with the greatest potential to be exploited is solar PV generation [40]. In this sense, in order to carry out this research, two Colombian regions, namely Medellín and Capurganá, were considered for the implementation of an EMS. These regions were selected for the following reasons:

i.

Each region has different climatic conditions. Medellín is a city located at latitude 6.2518${}^{°}$ N and longitude 75.5636${}^{°}$ W (Figure 1) and has an average annual temperature of 23 ${}^{°}$C. Likewise, the municipality of Capurganá is located at latitude 8.6167${}^{°}$ N and longitude 77.3333${}^{°}$ W (Figure 1) and has an average annual temperature of 28 ${}^{°}$C.

ii.

The selected regions have different PV power potentials. In Figure 1, it can be highlighted that the solar energy potential of Medellín is approximately 4.6 kWh/m${}^2$/day, while in Capurganá, the solar energy potential is approximately 3.7 kWh/m${}^2$/day [41].

iii.

Due to their locations, the energy consumption habits of the two Colombian regions are very different [42]. Medellín is a city that belongs to the SIN and therefore has access to electricity 24 h a day. Its economic activities are based on industry and commerce, and it is also the second most important and populated city in the country. On the other hand, the municipality of Capurganá is a ZNI located in a place of difficult access, so its electricity generation is based on diesel, with an average of 19 h of access to electricity per day. Thus, its economic activities are based on fishing and agriculture.

As shown in Figure 1, these two regions are completely different in terms of solar generation and energy consumption habits. In addition, due to the fact that the country is divided into two large regions (from the point of view of energy infrastructure) based on how easy (i.e., regions belonging to the SIN) or difficult (i.e., regions belonging to the ZNI) it is to supply the energy demand of end users, the authors of this document have proposed a mathematical model that considers the variations in power generation and demand in order to evaluate the impact of implementing an EMS on the technical, economic, and environmental indicators of distribution systems located in urban areas and rural areas (i.e., Medellín and Capurganá, respectively).

1.5. Document Structure

This research paper is structured as follows: Section 2 presents the mathematical formulation of the operation problem regarding PV generators in distribution grids; Section 3 presents the EMS designed through a master–slave methodology that integrates the ALO with the successive approximations power flow method; Section 4 shows the generation and demand curves used to implement the proposed optimization strategy; Section 5 presents the main characteristics of the 27- and 33-node test systems used for validation as well as the parametric information necessary to calculate the value of each fitness function; Section 6 shows the results obtained for the operation of PV generators in each proposed simulation scenario; and Section 7 presents the conclusions and future works derived from this research project.

2. Optimal Operation of PV Generators

The problem regarding the optimal operation of PV generators in distribution grids can be represented through a nonlinear programming (NLP) model, where the decision variables of the problem (i.e., variables of a continuous nature) are associated with the injection of power by the PV generators, while nonlinearities of the model appear in the power flow formulation, given the nonlinear and non-convex nature of its equations [43]. This section presents the formulation of the objective functions and the set of constraints that represent the problem of optimal PV generator operation in distribution grids.

2.1. Formulation of the Objective Function

In distribution grids with PV generators, the objective function may be aimed at improving technical, economic, or environmental conditions, through which the following may be sought: minimizing power or energy losses, system operating costs, or CO${}_2$ emissions; improving voltage profiles; or reducing the loadability of distribution lines [44]. The selection of said objective function only depends on the needs and requirements of the network operator.

The objective functions selected for this research work were the minimization of operating costs, energy losses, and CO${}_2$ emissions, aiming to bring the distribution grid to an optimal point of operation while ensuring that the demand is met as economically as possible and with a low environmental impact.

2.1.1. Economic Indicator

The objective function considered in this study to improve the economic indices of the distribution grid is the minimization of the operating costs of the network, which is composed of the costs associated with the generation or purchase of energy at the main supply node (i.e., slack node or substation node) and the operation and maintenance costs of the PV generators installed within the [44] system. Each component of the objective function is presented in Equations (1)–(3).

$$min{E}_{cost}=f_1+f2,$$

(1)

$$f_1=C_{kWh}\left(\sum _{h\in \mathcal{H}}\sum _{i\in \mathcal{N}}{p}_{i,h}^s\Delta h\right),$$

(2)

$$f_2=C_{O\&M}\left(\sum _{h\in \mathcal{H}}\sum _{i\in \mathcal{N}}{p}_{i,h}^{pv}\Delta h\right),$$

(3)

where $E_{cost}$ represents the total operating costs of the MG; $f_1$ is the component of the objective function related to the purchase of energy at the terminals of the substation node; $f_2$ is the component of the objective function associated with the operation and maintenance costs of the PV generators; $C_{kWh}$ refers to the average cost of purchasing power at the substation node; $p_{i,h}^s$ is the active power provided by a conventional generator connected to a node i during a period of time h; $\Delta h$ is the duration of the period of time in which the electrical variables are assumed to be constant; $C_{O\&M}$ is the maintenance and operation cost of a PV generator; $p_{i,h}^{pv}$ is the active power generated by a PV source connected to a node i during a period of time h; and $\mathcal{N}$ and $\mathcal{H}$ are the sets that contain all the network nodes and the time periods in a daily operation scenario, respectively.

2.1.2. Technical Indicator

To improve the technical indicators of the MG, the minimization of network energy losses is taken as the objective function. This function goal is presented in Equation (4) [44].

$$min{E}_{loss}=\sum _{h\in \mathcal{H}}\sum _{l\in \mathcal{L}}{R}_l{I}_l^2\Delta h,$$

(4)

where the value of $E_{loss}$ represents the energy losses of the MG; $R_l$ is the resistance associated with a line l; and $I_l$ is the magnitude of the current flowing through the section of the network l. Note that $\mathcal{L}$ is the set containing all the distribution lines of the MG.

2.1.3. Environmental Indicator

Finally, aiming to improve the environmental indicators of the MG, the minimization of CO${}_2$ emissions is taken as the objective function. This objective function comprises the emissions associated with power generation at the substation node and those related to power generation at the PV generators [45]. This objective function is shown in Equation (5).

$$min{E}_{C{O}_2}=C{E}_s\left(\sum _{h\in \mathcal{H}}\sum _{i\in \mathcal{N}}{p}_{i,h}^s\Delta h\right),$$

(5)

where $E_{C{O}_2}$ represents the total emissions due to the operation of the distribution grid and $C{E}_s$ refers to the CO${}_2$ emission factor associated with conventional generation sources.

2.2. Set of Constraints

The set of constraints corresponds to the different operational limitations that can be found in electrical systems within a distributed generation environment, such as the balance of active and reactive power in each node, the voltage regulation limits, and the thermal limit of the drivers present in the network, among others [46]. In this vein, the complete list of constraints is presented below. Note that the set of restrictions is presented as a function of the real and imaginary parts of the nodal voltage and the current that circulates through the network sections.

Equality constraints (6) and (7) represent the balance of active and reactive power for each node of the system in each period of time.

$$p_{i,h}^s+p_{i,h}^{pv}-P_{i,h}^d=\sum _{l\in \mathcal{L}}{A}_{i,l}\left(V_{i,h}^{re}{I}_{l,h}^{re}+V_{i,h}^{im}{I}_{l,h}^{im}\right),\left\{\begin{array}{c}\forall i\in \mathcal{N},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(6)

$$q_{i,h}^s-Q_{i,h}^d=-\sum _{l\in \mathcal{L}}{A}_{i,l}\left(V_{i,h}^{re}{I}_{l,h}^{im}-V_{i,h}^{im}{I}_{l,h}^{re}\right),\left\{\begin{array}{c}\forall i\in \mathcal{N},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(7)

where $P_{i,h}^d$ and $Q_{i,h}^d$ are the active and reactive power demanded at a node i in a period of time h, respectively; $q_{i,h}^s$ is the reactive power generated by each conventional generator connected to a node i in a period of time h; $V_{i,h}^{re}$ and $V_{i,h}^{im}$ are the real and imaginary parts of the voltage at a node i during a period of time h; $I_{i,h}^{re}$ and $I_{i,h}^{im}$ are the real and imaginary parts of the current that circulates through a distribution line l during a period of time h; and $A_{i,l}$ is the node-to-branch incidence matrix and can be constructed as shown below [47].

• $A_{il}=1$ if line l is connected to the node i and the current flow is leaving this node.
• $A_{il}=-1$ if line l is connected to the node i and the current flow arrives at this node.
• $A_{il}=0$ if line l is not connected to the node i.

On the other hand, inequality constraints (8) and (9) define the lower and upper limits for the injection of active and reactive power by conventional generators. Similarly, (10) is an inequality constraint that defines the lower and upper limits of active power generation for the PV generators present in the MG.

$$P_i^{s,min}\le p_{i,h}^s\le P_i^{s,max},\left\{\begin{array}{c}\forall i\in \mathcal{N},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(8)

$$Q_i^{s,min}\le q_{i,h}^s\le Q_i^{s,max},\left\{\begin{array}{c}\forall i\in \mathcal{N},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(9)

$$P_i^{pv,min}\le p_{i,h}^{pv}\le P_i^{pv,max},\left\{\begin{array}{c}\forall i\in \mathcal{N},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(10)

where $P_i^{s,min}$ and $P_i^{s,max}$ are the active power limits associated with each conventional generator connected to a node i; $Q_i^{s,min}$ and $Q_i^{s,max}$ are the reactive power limits associated with each conventional generator connected to a node i; and $P_i^{pv,min}$ and $P_i^{pv,max}$ are the active power limits associated to each PV generator connected to a node k.

Equality constraints (11) and (12) define the real and imaginary parts of the currents that circulate in the distribution lines.

$$I_{l,h}^{re}=\frac{1}{R_l^2+X_l^2}\sum _{i\in \mathcal{N}}{A}_{i,l}\left(R_l{V}_{i,h}^{re}+X_l{V}_{i,h}^{im}\right),\left\{\begin{array}{c}\forall l\in \mathcal{L},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(11)

$$I_{l,h}^{im}=\frac{1}{R_l^2+X_l^2}\sum _{i\in \mathcal{N}}{A}_{i,l}\left(R_l{V}_{i,h}^{im}-X_l{V}_{i,h}^{re}\right),\left\{\begin{array}{c}\forall l\in \mathcal{L},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(12)

where $X_l$ is the reactance value associated with the distribution line l.

Equality constraints (13) and (14) present the voltage magnitude at a node i and the current magnitude in a distribution line l, respectively.  

$$V_{i,h}=\sqrt{{\left(V_{i,h}^{re}\right)}^2+{\left(V_{i,h}^{im}\right)}^2},\left\{\begin{array}{c}\forall i\in \mathcal{N},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(13)

$$I_{l,h}=\sqrt{{\left(I_{l,h}^{re}\right)}^2+{\left(I_{l,h}^{im}\right)}^2},\left\{\begin{array}{c}\forall l\in \mathcal{L},\\ \forall h\in \mathcal{H}\end{array}\right\}.$$

(14)

Finally, box-type constraints (15) and (16) define the lower and upper limits of voltage regulation for all nodes as well as the current capacity of all the distribution lines.

$$V_i^{min}\le V_{i,h}\le V_i^{max},\left\{\begin{array}{c}\forall i\in \mathcal{N},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(15)

$$0\le I_{l,h}\le I_l^{max},\left\{\begin{array}{c}\forall l\in \mathcal{L},\\ \forall h\in \mathcal{H}\end{array}\right\}.$$

(16)

Note that constraint (16) ensures that the thermal limit of each conductor present in the distribution grid is respected. With the power flow formulation presented in this section, the telescopic characteristic of electrical networks can be considered: as the distribution lines move away from the substation bus, their caliber is reduced in addition to the maximum current that the conductors can withstand [48].

Traditionally, PV generators are operated with a constant and pre-set generation factor, which is preserved during the optimization process [49,50,51]. This means that the PV generators installed in a distribution grid supply all the power they can generate to the grid (i.e., maximum power point tracking). One of the main disadvantages of this solution methodology is that it can negatively impact the operating conditions of the network in scenarios where PV power generation exceeds power demand [52,53]. To avoid this, the solution methodology presented in this paper takes advantage of freeing the PV generators, so that they are not forced to inject all the power they can generate in each period of time (i.e., deactivating maximum power point tracking), thus making it possible for the PV generator to operate optimally and ensuring the balance between generation and demand in each period of time [54]. This behavior can be modeled mathematically, as shown in (17).

$$p_{i,h}^{pv}\le P_i^{pv}{C}_h^{pv},\left\{\begin{array}{c}\forall i\in \mathcal{N},\\ \forall h\in \mathcal{H}\end{array}\right\},$$

(17)

where $P_i^{pv}$ is the nominal power of the PV generator located at a node i and $C_h^{pv}$ is the expected PV generation behavior curve for the area where the distribution grid is located.

Inequality constraints (10) and (17) show that the output power of the PV generator in a period of time h can take values between the lower limit of power and the product between the nominal power and the solar generation curve.

Note that the NLP model defined from (1)–(17) is the general representation of the problem regarding the optimal operation of PV generators in distribution grids. However, due to the active and reactive power flow constraints, this problem is nonlinear and non-convex, so optimization techniques and numerical methods must be used to solve it [55]. Consequently, a master–slave methodology combining the antlion optimizer (ALO) and the successive approximations power flow method is proposed in this research.

3. Master–Slave Solution Methodology

To solve the operation problem regarding PV generators in distribution grids that was modeled in the previous section, this study proposes the implementation of an EMS designed with a master–slave methodology that uses the ALO as the master stage [56] and the successive approximations method as the slave stage [57]. The master stage defines the power that each PV generator must dispatch for each hour of operation, while the slave stage deals with the constraints associated with the power flow and defined from (6) to (17). Next, the coding used to represent the problem will be described as well as each of the components of the proposed methodology.

3.1. Proposed Coding

The structure of the coding to adapt to the problem under study is presented in Figure 2.

This figure shows a vector of size $1\times (N_{pv}·H)$, where each column of the vector corresponds to the power generated by each PV generator in a given period of time, where $N_{pv}$ is the number of PV generators available in the distribution grid and H is the hours of availability of the energy resource for each generator. In this vector, each PV generator has 13 power values to supply to the grid within the period of time in which the solar resource is available in Colombia [58]. Therefore, to represent this behavior, a vector of $13\times N_{pv}$ elements is required for $N_{pv}$ PV generators. In the example of Figure 2, PV generator 1 generates a power of 0.01 MW during time period 1, 0.12 MW during time period 2, 0.25 MW during time period 3, and 0.18 MW during time period 13. Similarly, the PV generator $N_{pv}$ generates a power of 0.15 MW during time period 1, 0.20 MW during time period 2, 0.35 MW during time period 3, and 0.05 MW during time period 13. Note that the technical limits of PV power generation must be established in advance to determine the minimum and maximum generation power range assigned to each of the time periods in which the PV generators will operate. In addition, note that this coding works adequately with regard to establishing the optimal operation of the existing PV generators in distribution grids, as it allows considering the hourly variation of generation based on the power demanded.

3.2. Master Stage: Antlion Optimization

The ALO is a bio-inspired metaheuristic optimization technique based on the hunting mechanism of antlions in nature [56,59]. The name of this insect originates from its hunting tactics and its favorite prey (i.e., ants). Its hunting mechanism is based on using its powerful jaw to dig a cone-shaped hole in the ground [60]. After digging a cone-shaped trap, these insects hide at the bottom of the pit and wait for their prey to be trapped. Once trapped, antlions will attack their prey until they eat it. After feeding, they clean and prepare the den for the next hunt. This behavior can be modeled mathematically through some simple rules of evolution, which will be explained below [56,59,60,61].

3.2.1. Initial Population

The initial populations of ants and antlions take the structure shown in (18) and (19), respectively.

$$A^t=\left[\begin{array}{cccc}{a}_{11}^t& a_{12}^t& \cdots & a_{1,N_v}^t\\ a_{21}^t& a_{22}^t& \cdots & a_{2,N_v}^t\\ ⋮& ⋮& \ddots & ⋮\\ a_{N_{i}1}^t& a_{N_{i}2}^t& \cdots & a_{N_i,N_v}^t\end{array}\right],$$

(18)

$$A{l}^t=\left[\begin{array}{cccc}a{l}_{11}^t& a{l}_{12}^t& \cdots & a{l}_{1,N_v}^t\\ a{l}_{21}^t& a{l}_{22}^t& \cdots & a{l}_{2,N_v}^t\\ ⋮& ⋮& \ddots & ⋮\\ a{l}_{N_{i}1}^t& a{l}_{N_{i}2}^t& \cdots & a{l}_{N_i,N_v}^t\end{array}\right],$$

(19)

where $A^t$ and $A{l}^t$ are the arrays that store the position of the ants and the position of the antlions’ hiding places in iteration t, respectively; $N_i$ is the number of individuals that make up each population; and $N_v$ is the number of variables or the dimension of the solution space, i.e., the product between the number of PV generators available in the distribution grid and the hours of availability of the solar resource ($N_{pv}·13$). To create an initial population of individuals that are capable of maintaining the structure shown in Figure 2, Equation (20) is used, which generates a matrix X of random numbers within the operating range of PV generators.

$$X=X^{min}ones(N_i,N_v)+(X^{max}-X^{min})rand(N_i,N_v),$$

(20)

where $ones(N_i,N_v)$ is a matrix filled by ones; $rand(N_i,N_v)$ is an array of random numbers that can take values between 0 and 1 and are generated by means of a uniform distribution; and $X^{min}$ and $X^{max}$ are vectors representing the lower and upper bounds of the solution space, as shown below:

$$X^{min}=\left[\begin{array}{ccc}{X}_1^{min}& ,\cdots ,& X_{N_{pv}}^{min}\end{array}\right],$$

$$X^{max}=\left[\begin{array}{ccc}{X}_1^{max}& ,\cdots ,& X_{N_{pv}}^{max}\end{array}\right],$$

with $X_i^{min}$ and $X_i^{max}$ being the vectors that contain the lower and upper limits of the decision variables associated with the dispatch of a PV generator i, as shown below:

$$X_i^{min}=\left[\begin{array}{ccc}{x}_{i,1}^{min}& ,\cdots ,& x_{i,H}^{min}\end{array}\right],$$

$$X_i^{max}=\left[\begin{array}{ccc}{x}_{i,1}^{max}& ,\cdots ,& x_{i,H}^{max}\end{array}\right].$$

Once the initial population of ants and antlions has been generated, the function of each of the individuals is evaluated, as shown in (21) and (22), respectively.

$$F_f\left(A^t\right)=\left[\begin{array}{c}{F}_f\left(A_1^t\right)\\ F_f\left(A_2^t\right)\\ ⋮\\ F_f\left(A_{N_i}^t\right)\end{array}\right]$$

(21)

$$F_f\left(A{l}^t\right)=\left[\begin{array}{c}{F}_f\left(A{l}_1^t\right)\\ F_f\left(A{l}_2^t\right)\\ ⋮\\ F_f\left(A{l}_{N_i}^t\right)\end{array}\right]$$

(22)

During this process, the best antlion is selected as the best solution found so far ($A{l}_{best}$).

Remark 1.

$F_f(·)$ represents the adaptation function to be minimized. In this research document, the minimization of operating costs (see (1)), the minimization of energy losses (see (4)), or the minimization of emissions of $C{O}_2$ (see (5)) can be selected according to the needs of the network operator.

3.2.2. Building the Trap

In the ALO, each ant can fall into only one trap in each iteration. It is for this reason that the roulette operator is used to model the hunting capacity of antlions (for more information, see [56]). By means of this operator, the antlion is selected according to its adaptation function value, which allows selecting the most suitable antlion for hunting [59].

3.2.3. The Ants Slide Toward the Antlion

Considering the above, it is possible to deduce that antlions are capable of building traps of a quality proportional to their adaptation function. In this sense, when an ant falls into a trap, the antlion starts shooting sand from the center of the pit to prevent the ant from escaping [59,61]. This behavior is modeled mathematically by reducing the space of solutions an ant can move through [56]. Therefore, the lower and upper bounds of all decision variables are reduced and updated, as shown in (23) and (24).

$$c^t=\frac{X^{min}}{R}$$

(23)

$$d^t=\frac{X^{max}}{R},$$

(24)

where $c^t$ and $d^t$ are the modified vectors of the lower and upper bounds of the solution space at iteration t, and R is the radius defined in (25).

$$R={10}^w\frac{t}{t_{max}},$$

(25)

where $t_{max}$ is the maximum number of iterations defined for the exploration and exploitation of the solution space, and w is a constant that takes a value according to the current iteration, as can be seen in (26).

$$w=\left\{\begin{array}{cc}\hfill 2\hfill & \hfill \mathrm{if}t>0.1t_{max}\hfill \\ \hfill 3\hfill & \hfill \mathrm{if}t>0.5t_{max}\hfill \\ \hfill 4\hfill & \hfill \mathrm{if}t>0.75t_{max}\hfill \\ \hfill 5\hfill & \hfill \mathrm{if}t>0.9t_{max}\hfill \\ \hfill 6\hfill & \hfill \mathrm{if}t>0.95t_{max}\hfill \end{array}\right.$$

(26)

The constant w is in charge of adjusting the level of precision of the exploitation in addition to ensuring the convergence of the algorithm. This is due to the fact that as the iterations increase, the space of the solutions is reduced [60].

3.2.4. Trapping the Ants in the Antlion Pits

As seen in the previous step, the traps limit the solution space through which the ants can move. To mathematically model this behavior, the lower and upper bounds of the solution space are adjusted in each iteration, so that the ant moves around the previously selected antlion trap, as shown in Equations (27) and (28).

$$c^t=\left\{\begin{array}{cc}\hfill A{l}_i^t+c^t\hfill & \hfill \mathrm{if}rand<0.5\hfill \\ \hfill A{l}_i^t-c^t\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$$

(27)

$$d^t=\left\{\begin{array}{cc}\hfill A{l}_i^t+d^t\hfill & \hfill \mathrm{if}rand<0.5\hfill \\ \hfill A{l}_i^t-d^t\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.,$$

(28)

where $A{l}_i^t$ is the antlion i selected by the roulette operator in iteration t, and $rand$ is a random number between 0 and 1.

3.2.5. Random Ant Walks

On the other hand, ants move randomly throughout the solution space in search for food. It is for this reason that random walks are used to represent their movement, which can be expressed as shown in (29) [56].  

$$\begin{array}{cc}\hfill RW=[& \hfill 0,cumsum(2r^1-1),cumsum(2r^2-1),\cdots ,\hfill \\ \hfill & \hfill cumsum(2r^t-1),\cdots ,cumsum(2r^{t_{max}}-1)],\hfill \end{array}$$

(29)

where $RW$ is the vector containing the positions of the ants’ random walk, $cumsum$ computes the cumulative sum, and $r^t$ is a stochastic function defined as:

$$r^t=\left\{\begin{array}{cc}\hfill 1\hfill & \hfill \mathrm{if}rand>0.5\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$$

(30)

According to Equation (29), the ants randomly update their position randomly. However, due to the effect of Equations (27) and (28), as the iterations progress, the space of solutions through which an ant can move is reduced. Therefore, (31) must be applied so that random walks occur within the solution space.

$$R{W}^t=\frac{\left(R{W}^t-a\right)\left(d^t-c^t\right)}{\left(b-a\right)}+c^t,$$

(31)

where a and b represent the minimum and maximum values of the random walk.

3.2.6. Elitism

This step is one of the most important steps in the ALO, as it allows keeping the best solutions obtained during the optimization process within the antlion population [56,59]. In this sense, the best solution obtained so far ($A{l}_{best}$) is capable of affecting the way in which the ants move within the solution space as the iterations advance. Therefore, the ant walks randomly around the roulette-selected antlion, which is the best antlion found so far, as shown in (32).

$$A_i^t=\frac{R{W}_A^t+R{W}_{best}^t}{2},$$

(32)

where $A_m^t$ represents the position of an ant m at iteration t. $R{W}_i^t$ is the random walk around the antlion selected by the roulette wheel at iteration t, while $R{W}_{best}^t$ is the random walk around the best antlion found so far in iteration t.

3.2.7. Catching Prey and Rebuilding the Pit

In the last stage of the antlion’s hunt, the ant has fallen to the bottom of the pit and has been caught in its jaw. Later, the antlion eats its prey and cleans and prepares the pit to hunt its next prey. The ALO replicates this process by assuming that the objective function value of the ants’ new position is better than that of the antlion selected by the roulette mechanism [56,59]. In this sense, the antlion updates its position, taking the position of the hunted ant to increase its chances of finding new prey. This behavior can be modeled as shown in (33).

$$A{l}_i^{t+1}=\left\{\begin{array}{cc}\hfill A_m^{t+1}\hfill & \hfill \mathrm{if}{F}_f\left(A_m^{t+1}\right)<F_f\left(A{l}_i^t\right)\hfill \\ \hfill A{l}_i^t\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$$

(33)

Similarly, the position of the best antlion found so far should be updated if there is an antlion whose fitness function is better than that of the former, as shown by the expression (34).

$$A{l}_{best}^{t+1}=\left\{\begin{array}{cc}\hfill A_i^{t+1}\hfill & \hfill \mathrm{if}{F}_f\left(A{l}_i^{t+1}\right)<F_f\left(A{l}_{best}^t\right)\hfill \\ \hfill A{l}_{best}^t\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$$

(34)

3.2.8. Stopping Criteria

It is common for metaheuristic algorithms to apply a stopping criterion in order to decide whether the algorithm stops or continues in the next iteration [59]. This is completed in order to reduce the computation time it takes for the algorithm to arrive at a solution. In this research paper, the following two stopping criteria for the ALO are proposed:

• Maximum number of iterations: the execution of the ALO will stop when the iteration counter t reaches a maximum number ($t_{max}$).
• Number of non-improvement iterations: the execution of the ALO will stop when the best antlion found so far is not updated after $t_{non}$ consecutive iterations.

3.3. Slave Stage: Successive Approximations Power Flow Method

By using the successive approximations method, it is possible to iteratively solve the active and reactive power balance equations described in (6) and (7), respectively. This numerical method allows the slave stage to determine the value of the objective function for each of the individuals that make up the population, evaluating the set of restrictions presented in the NLP model (Section 2). Considering the above, the following is the power flow formulation using successive approximations, as originally described by [57]. These authors take advantage of the characteristics of distribution systems (where there is a single $slack$ node and $n-1$ loading nodes) to arrive at the recursive formula shown in (35).

$${\mathbb{V}}_{d,h}^{t+1}=-{\mathbf{Y}}_{dd}^{-1}\left[{\mathbf{diag}}^{-1}\left({\mathbb{V}}_{d,h}^{t,*}\right)({\mathbb{S}}_{d,h}^{*}-{\mathbb{S}}_{pv,h}^{*})+{\mathbf{Y}}_{ds}{\mathbb{V}}_{s,h}\right],$$

(35)

where t is the iterative counter; ${\mathbb{V}}_d$ is the vector containing the voltages at the demand nodes; ${\mathbf{Y}}_{dd}$ is the component of the admission matrix that associates the demand nodes with each other; ${\mathbb{S}}_d$ corresponds to the complex power demanded at all load nodes; ${\mathbb{S}}_{pv,h}$ is the vector that contains the active power for each PV generator in each time period h; ${\mathbf{Y}}_{ds}$ is the component of the admission matrix that associates the demand nodes with the $slack$ node; and ${\mathbb{V}}_s$ is the vector containing the voltage at the $slack$ node. Note that the value of ${\mathbb{S}}_{pv}$ is provided by the master stage and is a vector that respects the encoding shown in Figure 2. Note that the iterative process ends when the maximum difference of the demand voltage magnitudes between two consecutive iterations is less than a maximum admissible error (i.e., convergence criterion), as shown in (36):

$$max\left\{\begin{array}{c}\left|\right|{\mathbb{V}}_{d,h}^{t+1}|-|{\mathbb{V}}_{d,h}^t\left|\right|\end{array}\right\}\le ϵ,$$

(36)

where $ϵ$ is the convergence error, which, for this research paper, has a value of $1\times {10}^{-10}$, as recommended by the authors of [62].

Once the power flow has been solved using the successive approximations method for all the time periods h, it is possible to calculate the value of $E_{loss}$. To this effect, it is necessary to calculate the current that circulates through the distribution lines in each period of time h, as shown in (37) [63].

$${\mathbb{I}}_{l,h}={\mathbf{Y}}_p{\mathbf{A}}^T\left[\begin{array}{c}{\mathbb{V}}_{s,h}\\ {\mathbb{V}}_{d,h}\end{array}\right],$$

(37)

where ${\mathbb{I}}_{l,h}$ is the vector in the complex domain that contains the current flowing through the distribution lines of the system; ${\mathbf{Y}}_p$ is the primitive admission matrix containing the inverse of the impedance of each line on its diagonal; and $\mathbf{A}$ is the incidence matrix. Likewise, to obtain the value of $E_{cost}$ and $E_{C{O}_2}$, it is necessary to calculate the power generated at the terminals of the conventional generator, as shown in (38).

$${\mathbb{S}}_{s,h}^{*}=\mathbf{diag}\left({\mathbb{V}}_{s,h}^{*}\right)\left({\mathbf{Y}}_{ss}{\mathbb{V}}_{s,h}+{\mathbf{Y}}_{sd}{\mathbb{V}}_{d,h}\right),$$

(38)

where ${\mathbb{S}}_{s,h}$ is the vector containing the complex power produced at node $slack$, and ${\mathbf{Y}}_{ss}$ is the component of the admission matrix associated with the $slack$ node.

In this way, it is possible to solve (37) and (38) and, consequently, to determine the value of the objective functions proposed in (1), (4) and (5). Additionally, in order to discard the potentially unfeasible solutions of the master stage, i.e., they violate the restrictions shown in (6)–(17), the objective functions described in (1), (4) and (5) are replaced by the adaptation functions shown in (39)–(41).  

$$F_{f1}=\left(\begin{array}{c}{E}_{cost}+{\alpha }_1\underset{h}{max}\left\{0,V_{i,h}-V_i^{max}\right\}\\ -{\alpha }_2\underset{h}{min}\left\{0,V_{i,h}-V_i^{min}\right\}\\ -{\alpha }_3\underset{h}{min}\left\{0,\mathrm{real}(p_{i,h}^s-P_i^{s,min})\right\}\\ +{\alpha }_4\underset{h}{max}\left\{0,I_{l,h}-I_l^{max}\right\}\end{array}\right),$$

(39)

$$F_{f2}=\left(\begin{array}{c}{E}_{loss}+{\alpha }_1\underset{h}{max}\left\{0,V_{i,h}-V_i^{max}\right\}\\ -{\alpha }_2\underset{h}{min}\left\{0,V_{i,h}-V_i^{min}\right\}\\ -{\alpha }_3\underset{h}{min}\left\{0,\mathrm{real}(p_{i,h}^s-P_i^{s,min})\right\}\\ +{\alpha }_4\underset{h}{max}\left\{0,I_{l,h}-I_l^{max}\right\}\end{array}\right),$$

(40)

$$F_{f3}=\left(\begin{array}{c}{E}_{C{O}_2}+{\alpha }_1\underset{h}{max}\left\{0,V_{i,h}-V_i^{max}\right\}\\ -{\alpha }_2\underset{h}{min}\left\{0,V_{i,h}-V_i^{min}\right\}\\ -{\alpha }_3\underset{h}{min}\left\{0,\mathrm{real}(p_{i,h}^s-P_i^{s,min})\right\}\\ +{\alpha }_4\underset{h}{max}\left\{0,I_{l,h}-I_l^{max}\right\}\end{array}\right),$$

(41)

where $F_{f1}$, $F_{f2}$, and $F_{f3}$ correspond to the value of the adaptation function related to the economic, technical, and environmental indicators, respectively; and ${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3$, and ${\alpha }_4$ are the penalty factors applied to each objective function. These penalty factors apply as long as the solutions provided by the master stage include the NLP model constraints described in Section 2. In this research paper, the value of the penalty factors is set as $1\times {10}^3$, where each factor has its corresponding units.

Algorithm 1 provides a general description of the procedure followed by the proposed master–slave methodology to solve the problem regarding the optimal operation of PV generators.

Algorithm 1: General implementation of the master–slave methodology to solve optimization problems.

4. Variable Power Demand and Generation

To implement the strategy, it is necessary to have the energy generation curves of the solar resource as well as the energy demand curves of the distribution system’s users. The main characteristic of the data that make up these curves is that they are variable and depend on weather conditions and consumption habits in the region where the network is located. This study considers the behavior of the user and the solar resource for two areas of Colombia in particular. The first zone is the city of Medellín, in the department of Antioquia, which is an urban zone [64]. The second area is the municipality of Capurganá, Chocó, which is a rural area [65].

4.1. Solar Generation Curves

The power generated at the terminals of the PV generation systems depends on the weather conditions associated with solar radiation and the temperature to which the solar panels are exposed, which makes this type of generation source non-dispatchable [66]. It is for this reason that in order to correctly operate electrical networks with the presence of PV generators in urban and rural areas, the power produced with the available solar resource must be accurately known. In the specialized literature, different formulations have been proposed to determine the power produced by the radiation and the temperature of the panel. The ones that stand out the most are those based on energy balance and on the efficiency of the solar panel [67]. However, according to [67], the most accurate formulations are those based on the energy balance of the panel, since the solar energy absorbed by a panel is converted into electrical energy and heat.

To calculate the output power of a PV generator, Equation (42) is used.

$$p_{i,h}^{pv}=P_i^{pv}{f}_{pv}\left(\frac{G_h^T}{G_i^{T,STC}}\right)\left[1+{\alpha }_p\left(T_{i,h}^c-T_i^{c,STC}\right)\right]$$

(42)

where $f_{pv}$ is a PV power reduction factor that considers external conditions that may affect the power production of a panel (i.e., impurities, shading, and reduced efficiency, among others); $G_h^T$ is the solar radiation that falls on a PV generator in a period of time h; $G_i^{T,STC}$ is the solar radiation of the PV generator located at node i under standard test conditions; ${\alpha }_p$ is the power temperature coefficient; $T_{i,h}^c$ is the surface temperature of the PV generator located at a node i during a period of time h; and $T_i^{c,STC}$ is the surface temperature of the PV generators located at a node i under standard testing conditions. Note that the surface temperature of the PV generator can be calculated as shown in (43).

$$T_{i,h}^c=T_h^a+G_h^T\left(\frac{T_i^{c,NOCT}-T_i^{a,NOCT}}{G_i^{T,NOCT}}\right)\left(1-\frac{{\eta }_i^c}{\tau \alpha }\right),$$

(43)

where $T_h^a$ is the ambient temperature to which the PV generator is exposed in a period of time h; $T_i^{c,NOCT}$ is the nominal surface temperature of the PV generator located at node i when exposed to radiation $G_i^{T,NOCT}$ at an ambient temperature of $T_i^{a,NOCT}$; ${\eta }_i^c$ is the electrical efficiency of the PV generator located at a node i; $\tau$ is the solar transmittance of the PV generator; and $\alpha$ is the solar absorption of the PV generator.

By analyzing Equations (42) and (43), it can be concluded that the output power of a PV generator is a function of solar radiation and the ambient temperature of the area where the distribution grid is located, since the rest of the variables are constant parameters that depend on the standard testing conditions and the type of material with which the panels of the PV generator are built.

To determine the behavior curve for power generation based on the solar resource available in the study areas (i.e., urban and rural), the parametric information shown in

Table 2. Parametric information of the PV generation sources.

Parameter	Value	Unit
$P_i^{pv}$	1	W
$f_{pv}$	0.95	-
$G_i^{T,STC}$	1000	W/m${}^2$
${\alpha }_p$	−0.0045	1/${}^{°}$C
$T_i^{c,STC}$	25	${}^{°}$C
$T_i^{c,NOCT}$	46	${}^{°}$C
$G_i^{T,NOCT}$	800	W/m${}^2$
$T_i^{a,NOCT}$	20	${}^{°}$C
${\eta }_i^c$	0.141	-
$\tau \alpha$	0.9	-

is used. This information was taken from [67,68], assuming that the material of the solar panels is polycrystalline silicon.

Remark 2.

Note that by setting the nominal power of the PV generators at 1 W, a generation curve is obtained which varies between 0 and 1. That is to say, a curve is obtained using a per-unit representation that denotes the behavior of the solar resource in the areas under study ($C_h^{pv}$), which serves as input to solve the NLP model presented in Section 2.

4.1.1. Urban Case: Medellín, Antioquia, Colombia

For the urban area, climatological data regarding solar radiation and ambient temperature were taken from the NASA database [69], which was created for predicting energy resources worldwide. Solar radiation and ambient temperature data for 2019 (i.e., from 1 January to 31 December) were taken with a 1 h sampling.

For this study, the collected solar radiation and ambient temperature data were averaged per hour, as shown in

Table 3. Solar radiation data (W/m${}^2$), ambient temperature (${}^{°}$C), and behavior (p.u.) for an average day in the regions under study.

Region	Medellín			Capurganá
Hour	$G_T$	$T_a$	$C_{\mathit{pv}}$	$G_T$	$T_a$	$C_{\mathit{pv}}$
1	0	16.14132	0	0	24.44252	0
2	0	15.90636	0	0	24.32474	0
3	0	15.68132	0	0	24.22545	0
4	0	15.46022	0	0	24.14674	0
5	0	15.27545	0	0	24.08422	0
6	0	15.10329	0	0	24.03482	0
7	46.02425	15.15718	0.04541	29.14570	24.10367	0.02770
8	190.83559	16.15636	0.18424	142.11066	24.78126	0.13277
9	362.83753	17.43868	0.34100	291.61926	25.68211	0.26622
10	526.64647	18.87312	0.48161	431.95384	26.63671	0.38547
11	640.99058	20.27438	0.57375	540.61581	27.47515	0.47362
12	709.05312	21.36342	0.62572	605.16362	28.10252	0.52397
13	701.86370	21.98721	0.61809	606.93027	28.46775	0.52442
14	626.82690	22.12107	0.55716	583.07479	28.56923	0.50519
15	499.86074	21.83071	0.45236	490.55904	28.42334	0.43065
16	346.26581	21.20351	0.32052	359.22033	28.03460	0.32148
17	186.66671	20.38668	0.17693	204.48775	27.44945	0.18722
18	52.33403	19.35951	0.05066	64.51775	26.69008	0.06034
19	0.50986	18.32258	0.00050	3.17460	25.89016	0.00300
20	0	17.72414	0	0	25.39227	0
21	0	17.29586	0	0	25.09285	0
22	0	16.96148	0	0	24.87663	0
23	0	16.67395	0	0	24.70841	0
24	0	16.40545	0	0	24.56926	0

. Thus, when applying (42) and (43) with the data in Table 3, together with the data presented in Table 2, the average behavior of the solar resource for a typical day in Medellín is obtained, as shown in Figure 3 and Table 3.

4.1.2. Rural Case: Capurganá, Chocó, Colombia

The same methodology was applied in the rural study area. The solar radiation and ambient temperature data reported in Table 3 were taken in order to obtain the average behavior of the solar resource for a typical day in Capurganá, as shown in Figure 3 and Table 3.

4.2. Demand Curves

Having the demand behavior of the users, it becomes of vital importance to solve the problem regarding the operation of PV generators in distribution grids, as it allows determining how the power injection level should behave at the substation node, the PV generators to meet the user demand, and the power losses associated with energy transportation. To identify the behavior of the user energy demand in urban and rural areas, the demand history reported by the network operators is used.

4.2.1. Urban Case: Medellín, Antioquia, Colombia

For the urban area, power consumption data were taken from the historical reports made by the network operator Empresas Públicas de Medellín (EPM) [70]. Consumption data from 1 January to 31 December 2019 were taken with a sampling of 1 h. As for the power generation curves, the data collected were averaged per hour as shown in

Table 4. Power consumption data (kW) and behavior (pu) for an average day in the regions under study.

Region	Medellín		Capurganá
Hora	$P_d$	$P_{d,\mathit{pu}}$	$P_d$	$P_{d,\mathit{pu}}$
1	1,012,876.20	0.65509	428.04117	0.84573
2	974,315.40	0.63015	409.76717	0.80962
3	951,768.01	0.61557	317.81654	0.62795
4	952,169.92	0.61583	256.70648	0.50720
5	996,601.97	0.64457	51.70864	0.10217
6	1,080,667.80	0.69894	11.05835	0.02185
7	1,135,234.91	0.73423	32.49553	0.06421
8	1,226,850.93	0.79348	62.77491	0.12403
9	1,303,895.33	0.84331	119.17381	0.23547
10	1,354,781.01	0.87622	281.26057	0.55572
11	1,417,860.03	0.91702	333.09429	0.65813
12	1,462,589.11	0.94595	358.36076	0.70805
13	1,459,381.62	0.94388	368.01140	0.72712
14	1,439,889.28	0.93127	369.70917	0.73048
15	1,430,823.70	0.92541	379.97901	0.75077
16	1,426,481.64	0.92260	388.65478	0.76791
17	1,404,019.24	0.90807	386.78365	0.76421
18	1,373,896.43	0.88859	395.19266	0.78083
19	1,463,002.74	0.94622	430.88177	0.85134
20	1,478,398.44	0.95618	464.61670	0.91800
21	1,415,579.31	0.91555	476.40313	0.94128
22	1,310,824.08	0.84779	473.67462	0.93589
23	1,187,930.28	0.76831	467.29281	0.92328
24	1,086,900.38	0.70297	452.18590	0.89344

. With the data consigned in this table, the average power consumption behavior for a typical day in Medellín was obtained, as shown in Figure 4 and Table 4.

4.2.2. Rural Case: Capurganá, Chocó, Colombia

Regarding the rural area, power consumption data were taken from the reports of historical events by the IPSE [71], which is in charge of monitoring and supervising the non-interconnected electrical areas in Colombia with the purpose of promoting, developing, and implementing solutions. As in the previous case, the data collected were averaged per hour, as shown in Table 4. Similarly, it was possible to obtain the average consumption behavior for a typical day in Capurganá, as shown in Figure 4 and Table 4.

Remark 3.

It is important to mention that this research was carried out with data on solar radiation, ambient temperature, and power consumption for 2019 because it was the last normal year for the electricity sector worldwide (i.e., before the pandemic caused by the virus SARS-CoV-2).

Figure 5 summarizes the general methodology for obtaining the average behavior curves for solar generation and the demand by the end users of the two regions under study. This can be completed regardless of the region as long as there is access to historical data on solar radiation, ambient temperature, and power demand.

4.3. Other Considerations for the Proposed Ems

To determine the value of the objective functions defined in (1) and (5), the parametric information shown in

Table 5. Parameters used to calculate the economic and environmental indicators.

Parameter	Value	Unit	Parameter	Value	Unit
$C_{kWh}^{Urban}$	0.1302	USD/kWh	$C{E}_s^{Urban}$	0.1644	kg/kWh
$C_{kWh}^{Rural}$	0.2913	USD/kWh	$C{E}_s^{Rural}$	0.2671	kg/kWh
$C_{O\&M}^{pv}$	0.0019	USD/kWh	-	-	-

is used.

Table 5 shows the costs of generating energy in urban and rural areas, as well as the costs associated with the maintenance of PV generation systems. The energy generation costs for the study areas were taken from the reports made by the network operators to the Unified Information System (SUI) in 2019 [72,73]. The operation and maintenance costs of the PV generators were taken from [74]. Similarly, the emission factors associated with power generation in urban and rural areas are also shown in this table. The emissions factor for the urban area is the factor established by XM for the interconnected electrical system, to which EPM [75] belongs. The emissions factor for the rural area is the one associated with diesel fuel and was taken from the database of the Emission Factors of Colombian Fuels (FECOC) [76].

Finally, the voltage regulation limits for electrical systems with a voltage level greater than 1 kV and less than 62 kV were defined as $+5$ and $-10\%$ of the nominal voltage, as established by Colombian Technical Standard (NTC) No. 1340 [77].

5. Test Systems

This section presents the main characteristics of the test systems used in order to validate the proposed methodology for the operation of PV generators in distribution grids in both urban and rural areas. The 33-node test system was selected as the urban power grid, while the 27-node test system was selected for the rural study area.

5.1. Urban Test Feeder: 33-Node System

The 33-node test system is a radial distribution network originally proposed by the authors of [78]. This system is composed of 32 distribution lines and has a single conventional generator located at node 1 (the substation bus), which operates at a nominal voltage of 12.66 kV. The electrical configuration of this test system is depicted in Figure 6, and

Table 6. Parametric information of the 33-node test system.

Line l	Node i	Node j	${\mathit{R}}_{\mathit{ij}}\left(\mathbf{\Omega }\right)$	${\mathit{X}}_{\mathit{ij}}\left(\mathbf{\Omega }\right)$	${\mathit{P}}_{\mathit{j}}$(kW)	${\mathit{Q}}_{\mathit{j}}$(kvar)	${\mathit{I}}_{\mathit{l}}^{max}$(A)
1	1	2	0.0922	0.0477	100	60	385
2	2	3	0.4930	0.2511	90	40	355
3	3	4	0.3660	0.1864	120	80	240
4	4	5	0.3811	0.1941	60	30	240
5	5	6	0.8190	0.7070	60	20	240
6	6	7	0.1872	0.6188	200	100	110
7	7	8	1.7114	1.2351	200	100	85
8	8	9	1.0300	0.7400	60	20	70
9	9	10	1.0400	0.7400	60	20	70
10	10	11	0.1966	0.0650	45	30	55
11	11	12	0.3744	0.1238	60	35	55
12	12	13	1.4680	1.1550	60	35	55
13	13	14	0.5416	0.7129	120	80	40
14	14	15	0.5910	0.5260	60	10	25
15	15	16	0.7463	0.5450	60	20	20
16	16	17	1.2890	1.7210	60	20	20
17	17	18	0.7320	0.5740	90	40	20
18	2	19	0.1640	0.1565	90	40	40
19	19	20	1.5042	1.3554	90	40	25
20	20	21	0.4095	0.4784	90	40	20
21	21	22	0.7089	0.9373	90	40	20
22	3	23	0.4512	0.3083	90	50	85
23	23	24	0.8980	0.7091	420	200	85
24	24	25	0.8960	0.7011	420	200	40
25	6	26	0.2030	0.1034	60	25	125
26	26	27	0.2842	0.1447	60	25	110
27	27	28	1.0590	0.9337	60	20	110
28	28	29	0.8042	0.7006	120	70	110
29	29	30	0.5075	0.2585	200	600	95
30	30	31	0.9744	0.9630	150	70	55
31	31	32	0.3105	0.3619	210	100	30
32	32	33	0.3410	0.5302	60	40	20

shows the parametric information of the system [78]. Similarly, this table shows the maximum current that can be supported by each section of the network in order to correctly evaluate the constraint stipulated in (16). These values are determined from the maximum currents calculated using the power flow for the time with maximum power demand and the NTC 2050, assuming that the conductors operate at a temperature of about 60 ${}^{°}$C.

It is important to mention that this system has three PV generators, all with a nominal power of 2400 kW, located at nodes 12, 15 and 31.

5.2. Rural Test Feeder: 27-Node System

The 27-node test system also has a radial structure and was originally reported in [79]. Likewise, this system is made up of 26 distribution lines and has a single conventional generator located at node 1 (the substation bus), which operates at a nominal voltage of 23 kV. The electrical configuration of this test system is illustrated in Figure 7, and

Table 7. Parametric information of the 27-node test system.

Line l	Node i	Node j	${\mathit{R}}_{\mathit{ij}}\left(\mathbf{\Omega }\right)$	${\mathit{X}}_{\mathit{ij}}\left(\mathbf{\Omega }\right)$	${\mathit{P}}_{\mathit{j}}$(kW)	${\mathit{Q}}_{\mathit{j}}$(kvar)	${\mathit{I}}_{\mathit{l}}^{max}$(A)
1	1	2	0.0140	0.6051	0	0	240
2	2	3	0.7463	1.0783	0	0	165
3	3	4	0.4052	0.5855	297.50	184.37	95
4	4	5	1.1524	1.6650	0	0	85
5	5	6	0.5261	0.7601	255.00	158.03	70
6	6	7	0.7127	1.0296	0	0	55
7	7	8	1.6628	2.4024	212.50	131.70	55
8	8	9	5.3434	3.1320	0	0	20
9	9	10	2.1522	1.2615	266.05	164.88	20
10	2	11	0.4052	0.5855	85.00	52.68	70
11	11	12	1.1524	1.6650	340	210.71	70
12	12	13	0.5261	0.7601	297.50	184.37	55
13	13	14	1.2358	1.1332	191.25	118.53	30
14	14	15	2.8835	2.6440	106.25	65.85	20
15	15	16	5.3434	3.1320	255.00	158.03	20
16	3	17	1.2942	1.1867	255.00	158.03	70
17	17	18	0.7027	0.6443	127.50	79.02	55
18	18	19	3.3234	1.9480	297.50	184.37	40
19	19	20	1.5172	0.8893	340	210.71	25
20	20	21	0.7127	1.0296	85.00	52.68	20
21	4	22	8.2528	2.9911	106.25	65.85	20
22	5	23	9.1961	3.3330	55.25	34.24	20
23	6	24	0.7463	1.0783	69.70	43.20	20
24	8	25	2.0112	0.7289	255.00	158.03	20
25	8	26	3.3234	1.9480	63.75	39.51	20
26	26	27	0.5261	0.7601	170	105.36	20

shows the parametric information of the system [79,80] and the current thermal limit. Note that this system has three PV generators, all with a nominal power of 2400 kW, located at nodes 5, 9 and 19.

6. Numerical Results and Discussion

This section presents the numerical results obtained with the proposed methodology. Here, the ALO is compared against the following methodologies, which have been selected due to their excellent results when it comes to solving the problem of optimal power flow in electrical distribution systems [81,82,83]: (i) the particle swarm optimization (PSO) approach [84], (ii) the Chu and Beasley genetic algorithm (CBGA) [85], and (iii) the vortex search algorithm (VSA) [86]. The NLP model that represents the addressed problem has been implemented and solved in MATLAB version 2022a using our own scripts on a Dell Precision 3450 workstation with an Intel(R) Core(TM) i9-11900 [email protected] and 64.0 GB RAM running Windows 10 Pro 64-bit.

Additionally, to validate the EMS proposed in this document, the following three simulation scenarios are proposed: operating the PV generators with a focus on minimizing the economic indicator ($S_1$), the technical indicator ($S_2$), and the environmental indicator ($S_3$).

6.1. Optimization of the Algorithms’ Parameters

In order to ensure a fair comparison between the ALO and the optimization algorithms selected for comparison, each of them was tuned to guarantee the best performance when solving the studied problem. For instance, to select the parameters, the CBGA was used with an initial population of 40 individuals and a maximum number of iterations of 400.

Table 8. Parameters of the optimization algorithms used to solve the problem of PV generator operation.

Algorithm	Parameter	Value	Range
ALO	Number of individuals ($N_i$)	162	[1–200]
	Number of iterations ($t_{max}$)	1048	[1–2000]
	Number of non-improvement iterations ($t_{non}$)	546	[1–1000]
cPSO	Number of individuals ($N_i$)	180	[1–200]
	Number of iterations ($t_{max}$)	1559	[1–2000]
	Number of non-improvement iterations ($t_{non}$)	417	[1–1000]
	Cognitive component ($C_1$)	1.1773	[0–2]
	Social component ($C_2$)	1.5640	[0–2]
	Maximum inertia ($W_{max}$)	0.5549	[0–1]
	Minimum inertia ($W_{min}$)	0.4377	[0–1]
CBGA	Number of individuals ($N_i$)	40	[1–200]
	Number of iterations ($t_{max}$)	1561	[1–2000]
	Number of non-improvement iterations ($t_{non}$)	1561	[1–1000]
	Number of random mutations ($nM$)	2	[0–$N_v$]
VSA	Number of individuals ($N_i$)	126	[1–200]
	Number of iterations ($t_{max}$)	1591	[1–2000]
	Number of non-improvement iterations ($t_{non}$)	312	[1–1000]
	Interval of radius reduction (x)	0.0655	[0–0.1]

presents a summary of the parameters tuned for each optimization algorithm as well as their value and the interval in which they were selected.

Likewise, the selected algorithms were executed 100 consecutive times to find the best, the average, and the worst value of the proposed adaptation functions. In the same way, this study calculated the standard percentage deviation and the average computation time that each algorithm takes to find the power vector which each existing PV generator must inject in both test systems.

6.2. Urban Zone Simulations

The results obtained after implementing the optimization methodology in the urban test system are shown in

Table 9. Numerical results in the 33-bus system for the urban area.

Scenario 1: Economical Index
Method	${\mathit{E}}_{\mathit{cost}}$(USD)	Reduction (%)	Time (s)	SDT (%)	${\mathit{V}}_{min}$(pu)/Node/Hour	${\mathit{V}}_{max}$(pu)/Node/Hour	${\mathit{I}}_{max}$(%)/Line/Hour
Bench. Case	9931.66	-	-	-	0.9084/18/20	1/1/All	94.4924/14/20
CBGA	7409.25	25.40	2.6258	0.4816	0.9084/18/20	1/1/All	97.7687/11/14
PSO	7317.89	26.32	31.4035	1.4559	0.9084/18/20	1/1/All	100/11/16
VSA	7276.05	26.74	42.4710	0.5786	0.9084/18/20	1/1/All	100/30/14
ALO	7220.09	27.30	141.9966	0.0088	0.9084/18/20	1/1/All	100/11/10
Scenario 2: Technical Index
Method	${\mathit{E}}_{\mathit{loss}}$(kWh)	Reduction (%)	Time (s)	SDT (%)	${\mathit{V}}_{min}$(pu)/Node/Hour	${\mathit{V}}_{max}$(pu)/Node/Hour	${\mathit{I}}_{max}$(%)/Line/Hour
Bench. Case	3379.07	-	-	-	0.9084/18/20	1/1/All	94.4924/14/20
CBGA	2346.00	30.57	2.6054	0.1918	0.9084/18/20	1/1/All	99.1658/14/14
PSO	2332.05	30.98	31.9980	0.1045	0.9084/18/20	1/1/All	94.4924/14/20
VSA	2331.61	30.99	37.9801	0.0034	0.9084/18/20	1/1/All	94.4924/14/20
ALO	2331.51	31.00	140.1495	0.0009	0.9084/18/20	1/1/All	94.4924/14/20
Scenario 3: Environmental Index
Method	${\mathit{E}}_{{\mathit{CO}}_{\mathbf{2}}}$(kg CO${}_{\mathbf{2}}$)	Reduction (%)	Time (s)	SDT (%)	${\mathit{V}}_{min}$(pu)/Node/Hour	${\mathit{V}}_{max}$(pu)/Node/Hour	${\mathit{I}}_{max}$(%)/Line/Hour
Bench. Case	12,541.22	-	-	-	0.9084/18/20	1/1/All	94.4924/14/20
CBGA	9309.57	25.77	2.6236	0.4513	0.9084/18/20	1/1/All	99.2224/11/10
PSO	9198.27	26.66	30.0663	1.2151	0.9084/18/20	1/1/All	100/30/14
VSA	9152.05	27.02	43.9550	0.3082	0.9084/18/20	1/1/All	100/11/14
ALO	9068.94	27.69	140.5316	0.0075	0.9084/18/20	1/1/All	100/11/12

. The information in this table is presented from left to right as follows: the methodology used, the value obtained for the evaluated function, the reduction obtained for each algorithm with respect to the base case, the computation time, the percentage standard deviation, the minimum voltage and the node and the time at which it occurs, the maximum voltage and the node and the time at which it occurs, and, finally, the maximum chargeability, its distribution line, and the time when it occurs. In the same way, Figure 8 shows the power injected by each PV generator, as obtained by the ALO, during daily operation in the urban area for each of the simulation scenarios evaluated.

The results obtained for the urban test system during one day of operation reveal the following:

✓

The proposed methodology achieves the best results regarding the evaluation of the adaptation function. In $S_1$, it reaches a response of 7220.09 USD, evidencing an improvement of 2711.57 USD with respect to the base case, 189.16 USD with respect to the CBGA, 97.81 USD with respect to the PSO, and 55.96 USD with respect to the VSA. In $S_2$, the ALO obtains a response of 2331.51 kWh, showing an improvement of 1047.56 kWh with respect to the base case, 14.50 kWh with respect to the CBGA, 0.54 kWh with respect to the PSO, and 0.10 kWh with respect to the VSA. Finally, in $S_3$, the proposed methodology achieves a response of 9068.94 kg of CO${}_2$, evidencing an improvement of 3472.28 kg of CO${}_2$ with respect to the base case, 240.63 kg of CO${}_2$ with respect to the CBGA, 129.33 kg of CO ${}_2$ compared to PSO, and 83.11 kg of CO${}_2$ when compared to the VSA.

✓

In $S_1$, all the methodologies used allow for reductions of more than 25% in comparison with the base case, with ALO reporting the highest value (27.30%). When this methodology is compared with the other metaheuristic algorithms, there is a reduction in operating costs of approximately 1.90% with respect to the CBGA, 0.98% with respect to PSO, and 0.56% with respect to the VSA. In $S_2$, all techniques allow obtaining reductions of more than 30% with respect to the base case, with ALO reporting the highest value (31%). The ALO achieves a reduction in energy losses of approximately 0.43% with respect to the CBGA, 0.02% with respect to PSO, and 0.01% with respect to the VSA. Moreover, in $S_3$, the evaluated methodologies allow for reductions of more than 25% with respect to the base case, with the proposed methodology showing the highest value, with 27.69%. When this methodology is compared to the other metaheuristic algorithms, there is a reduction in CO${}_2$ emissions of approximately 1.92% with respect to the CBGA, 1.03% when compared to the PSO, and 0.67% when compared to the VSA.

✓

Regarding the computation times, the CBGA, PSO, and the VSA are faster than the proposed ALO in the three simulation scenarios. The ALO takes approximately 141.9966 s in $S_1$, 140.1495 s in $S_2$, and 140.5316 s in $S_3$ to solve the PV generator operation problem. This shows that in order to solve a multidimensional (39-dimensional) NLP model with continuous variables (i.e., a solution space with infinite combinations), the proposed methodology takes less than 2.5 min to converge to an optimal solution. This allows grid operators in an urban area to implement an EMS that is capable of evaluating infinite combinations of PV power injection for one day of operation with low processing times in order to find the best solution from an economic, technical, and environmental point of view.

✓

As for the standard deviation, the superiority of the proposed ALO can be appreciated in all simulation scenarios. In $S_1$, it achieves a reduction of 5386.42% with respect to the CBGA, 16,485.55% with respect to PSO, and 6491.21% with respect to the VSA. In $S_2$, the proposed methodology obtains a reduction of 21,768.05% with respect to the CBGA, 11,808.53% with respect to PSO, and 282.49% with respect to the VSA. Finally, in $S_3$, the ALO reports a reduction of 5887.46% when compared to the CBGA, 16,019.70% when compared to the PSO, and 3988.20% when compared to the VSA.

✓

Likewise, it can be noted that for the three simulation scenarios, the results obtained with the ALO satisfactorily comply with voltage regulation, staying within $-10$ and $+5\%$ of the nominal voltage in all periods under analysis. The minimum voltage is found at node 18 and hour 20 (i.e., the time of maximum power demand in the urban area), with a value of 0.9084 p.u. Meanwhile, the maximum voltage is found at the slack node, with a value of 1 p.u. in all time periods. The proposed EMS allows managing and taking advantage of the solar resource, making it possible to inject PV power in time periods 7 to 19 without exceeding the minimum and maximum voltages of the distribution grid when PV power is not generated and there is a peak in power demand.

✓

Finally, the results obtained by the ALO show that in $S_1$, the maximum chargeability is reported in distribution line 11 at hour 10, with a value of 100%. In $S_2$, the maximum chargeability takes place in distribution line 14 at hour 20, with a value of 94.4924%. Moreover, in $S_3$, the maximum chargeability occurs in distribution line 11 at hour 12, with a value of 100%. The proposed optimization methodology allows for a smart operation of the PV generators, making it possible to respect the maximum bearable current for each conductor in the system while also allowing the efficient management of the solar resource between time periods 7 and 19.

The above demonstrates the efficiency and robustness of the ALO in solving the problem regarding the operation of PV generators in distribution systems with the aim of optimizing the system from an economic, technical, or environmental point of view. The methodology obtains the best performance; it has the best response and repeatability as well as low processing times (less than 2.5 min). This makes it the best option to address this problem in the urban test system, as it obtains an optimal solution for each simulation scenario which respects the technical-operating conditions of the network (i.e., minimum voltage, maximum voltage, and maximum withstand current per conductor).

Figure 8 shows the response obtained by the ALO for the operation of each PV generator in the three proposed simulation scenarios. This figure illustrates the maximum power that can be injected by each PV generator (black curve) according to the climate conditions presented in Section 4 for the city of Medellín, Antioquia, Colombia. Similarly, the power injected by each PV generator is illustrated (red, green, and blue curves).

Figure 8 also shows the injection of PV power for the three simulation scenarios. In time periods 7, 8, 17, 18 and 19, the PV generators inject the maximum possible power, as it is then that less than 20% of the solar resource is available (Figure 3). Similarly, the maximum possible power is injected because the demand is between 70 and 90% (Figure 4). In time periods 9 to 16, when the power demand is between 80 and 100%, the PV generators operate intelligently, injecting enough power to guarantee the power balance of the MG, as well as to ensure that the technical-operating conditions of the network are not violated, since, if the maximum possible power was generated (i.e., maximum power point tracking), the demanded power would be exceeded in periods 10 to 14 (up to 28% at hour 12).

Additionally, the injected power is greater in scenarios $S_1$ and $S_3$, with approximately 2183.7 kW being generated in the period of time with the highest solar radiation (i.e., hour 12), with which it is possible to cover 62% of the energy demand. In this case, the aim is to reduce the dependence of the system on the slack generator, which means that the network seeks to generate the least amount of power with the conventional generator and supply the greatest amount of demand through PV generators. With this, it is possible to greatly reduce the energy purchasing costs by the network operator and the CO${}_2$ emissions caused by the conventional generator. On the other hand, in $S_2$, approximately 1869.6 kW are generated during hour 12, which is enough to supply 53% of the demand. In this scenario, unlike the two previous cases, the aim is to inject the necessary power to improve the voltage profiles of the distribution grid, thus reducing energy losses to the maximum in one day of operation.

6.3. Rural Zone Simulations

The numerical results obtained after implementing the methodology in the rural test system are shown in

Table 10. Numerical results in the 27-bus system for the rural area.

Scenario 1: Economical Index
Method	${\mathit{E}}_{\mathit{cost}}$(USD)	Reduction (%)	Time (s)	SDT (%)	${\mathit{V}}_{min}$(pu)/Node/Hour	${\mathit{V}}_{max}$(pu)/Node/Hour	${\mathit{I}}_{max}$(%)/Line/Hour
Bench. Case	18,543.84	-	-	-	0.9664/10/21	1/1/All	91.4469/5/21
CBGA	12,282.02	33.77	1.8066	0.4548	0.9664/10/21	1.0013/9/9	99.2995/8/14
PSO	12,104.61	34.72	21.8934	0.6662	0.9664/10/21	1.0020/9/8	100/8/16
VSA	12,052.94	35.00	30.7303	0.2148	0.9664/10/21	1.0025/9/9	100/8/16
ALO	12,022.40	35.16	125.9021	0.0007	0.9664/10/21	1.0009/9/8	99.9987/8/15
Scenario 2: Technical Index
Method	${\mathit{E}}_{\mathit{loss}}$(kWh)	Reduction (%)	Time (s)	SDT (%)	${\mathit{V}}_{min}$(pu)/Node/Hour	${\mathit{V}}_{max}$(pu)/Node/Hour	${\mathit{I}}_{max}$(%)/Line/Hour
Bench. Case	691.15	-	-	-	0.9664/10/21	1/1/All	91.4469/5/21
CBGA	559.51	19.05	1.8311	0.0649	0.9664/10/21	1/1/All	91.4469/5/21
PSO	558.28	19.22	27.8036	0.0717	0.9664/10/21	1/1/All	91.4469/5/21
VSA	558.22	19.23	30.2881	0.0170	0.9664/10/21	1/1/All	91.4469/5/21
ALO	558.20	19.24	126.5206	0.0005	0.9664/10/21	1/1/All	91.4469/5/21
Scenario 3: Environmental Index
Method	${\mathit{E}}_{{\mathit{CO}}_{\mathbf{2}}}$(kg CO${}_{\mathbf{2}}$)	Reduction (%)	Time (s)	SDT (%)	${\mathit{V}}_{min}$(pu)/Node/Hour	${\mathit{V}}_{max}$(pu)/Node/Hour	${\mathit{I}}_{max}$(%)/Line/Hour
Bench. Case	17,005.21	-	-	-	0.9664/10/21	1/1/All	91.4469/5/21
CBGA	11,192.67	34.18	1.8213	0.5073	0.9664/10/21	1.0026/9/8	96.4822/8/10
PSO	11,064.72	34.93	21.0702	1.3339	0.9664/10/21	1.0009/9/8	99.9999/8/15
VSA	11,023.51	35.18	30.4519	0.1853	0.9664/10/21	1.0012/9/9	100/8/15
ALO	10,985.75	35.38	131.3081	0.0004	0.9664/10/21	1.0012/9/9	99.9995/8/15

, which presents the same parameters as Table 9. In the same way, as in the previous case study, Figure 9 shows the power generated by each PV generator, as obtained by the ALO during a daily operation in the rural area for each of the simulation scenarios.

The results shown in Table 10 reveal the following:

✓

In the rural test system, the proposed methodology achieves the best results regarding the evaluation of the adaptation function. In $S_1$, it reaches a response of 12,022.40 USD, evidencing an improvement of 6521.44 USD with respect to the base case, 259.62 USD with respect to the CBGA, 82.21 USD with respect to PSO, and 30.54 USD with respect to the VSA. In $S_2$, the ALO obtains a response of 558.20 kWh, i.e., an improvement of 132.95 kWh with respect to the base case, 1.30 kWh with respect to the CBGA, 0.080 kWh with respect to PSO, and 0.013 kWh with respect to the VSA. Finally, in $S_3$, it achieves a response of 10,0985.75 kg of CO${}_2$, which represents an improvement of 6019.45 kg of CO${}_2$ with respect to the base case, 206.91 kg of CO${}_2$ with respect to the CBGA, 78.97 kg of CO${}_2$ with respect to PSO, and 37.76 kg of CO${}_2$ when compared to the VSA.

✓

In $S_1$, all the methodologies allow for reductions of more than 33% with respect to the base case, with ALO being the methodology that reports the highest value, with 35.16%. There is a reduction in operating costs of approximately 1.40% when compared to the CBGA, 0.44% with respect to PSO, and 0.16% with respect to the VSA. In $S_2$, the studied methodologies allow obtaining a reduction of more than 18.5% with respect to the base case, with ALO reporting the highest value (19.24%). By comparing ALO with the other selected optimization algorithms, a reduction in energy losses of approximately 0.19% is obtained with respect to the CBGA, 0.012% with respect to PSO, and 0.002% with respect to the VSA. Finally, in $S_3$, there are reductions of more than 34% with respect to the base case, with the proposed methodology reporting the highest reduction (35.38%). This constitutes a reduction in CO${}_2$ emissions of approximately 1.22% when compared to the CBGA, 0.46% when compared to PSO, and 0.22% in comparison with the VSA.

✓

As for the computation time, the CBGA, PSO, and the VSA are faster than the proposed methodology in the three simulation scenarios. The ALO took an average time of 125.9021 s in $S_1$, 126.52 s in $S_2$, and 131.3081 s in $S_3$ to solve the PV generator operation problem. Similarly, it took less than 2.5 min to reach an optimal solution to a complex problem from the dimensional and solution space point of view.

✓

Regarding the standard deviation, it can be seen that in $S_1$, the ALO obtained a reduction of 68,355.35% when compared to the CBGA, 100,180.36% when compared to PSO, and 32,227.72% when compared to the VSA. In $S_2$, it achieved a reduction of 14,199.43% with respect to the CBGA, 15,699.79% with respect to PSO, and 3646.82% with respect to the VSA. Finally, in $S_3$, the ALO reported a reduction of 122,168.07% with respect to the CBGA, 321,417.59% with respect to the PSO, and 44,565.02% with respect to the VSA.

✓

Similarly, for the three simulation scenarios, the results obtained by the proposed methodology comply satisfactorily with voltage regulation. The minimum voltage was found for node 10 during hour 21 (i.e., the time of maximum power demanded in the rural area), with a value of 0.9664 p.u. Meanwhile, the maximum voltage was found for node 9 during hour 8, with a value of 1.0009 p.u.

✓

In $S_1$, the maximum chargeability was reported for distribution line 8 during hour 15, with a value of 99.9987%. In $S_2$, the maximum chargeability was observed in distribution line 5 during hour 21, with a value of 91.4469%. Finally, in $S_3$ the maximum chargeability ttook place in distribution line 8 during hour 15, with a value of 99.9995%.

The results indicate that the proposed ALO shows the best performance, namely a better response and repeatability, with low processing times (less than 2.5 min). This makes this solution methodology the best option to address the problem of operating PV generators in the rural test system, as it obtains an optimal solution for each simulation scenario which respects the technical-operating conditions of the network.

Finally, Figure 9 shows the response obtained by the ALO for the operation of each PV generator in the three simulation scenarios proposed for the rural test system. This figure illustrates the same information as in Figure 8. However, the maximum power that can be generated varies according to the weather conditions outlined in Section 4 for Capurganá, Chocó, Colombia.

Figure 9 shows the PV power injection for the three simulation scenarios. In time periods 17, 18, and 19, the PV generators inject the maximum possible power because less than 20% is available from the solar resource (Figure 3) and it is necessary to supply a power demand that is above 70% (Figure 4).

In $S_1$ and $S_3$, from time period 9, when the demand begins to progressively increase, the generators PV1 and PV3 try to inject a power equal to or very close to the maximum power that can be supplied in rural areas, while PV2 does not inject its maximum power; it generates enough power to ensure the power balance of the network and its technical-operating conditions. The three generators generate approximately 3039.8 kW during hour 14, covering 100% of the demand together with the power losses. Similarly, it can be noted that not only is the entire demand covered at hour 14, but power consumption and system losses are also covered from hours 8 to 13, thus allowing for significant operating costs and CO${}_2$ emissions reductions with respect to the base case.

In $S_2$, in time periods 7 to 16, when the power demand is between 0 and 80%, the PV generators manage the solar resource in such a way that they only provide enough power to guarantee the power balance of the distribution grid and ensure that the technical-operating conditions of the network are not violated. This evidences a progressive increase in PV generation as power demand increases, injecting approximately 1901.49 kW at hour 15 and managing to cover 62% of the power demanded by users. As in the urban test system, the operation of the PV generators seeks to inject the necessary power to improve the voltage profiles of the distribution grid and thus reduce the energy losses to the maximum for one day of operation.

6.4. Complementary Analysis and Discussion

Based on the results obtained in Section 6, the following can be highlighted: (i) by comparing the reductions in the operating costs of both test systems, a greater reduction is made evident in the rural area, i.e., 7.86% more than in the urban area. This represents savings of approximately 3809.87 USD in one day of operation for the Capurganá network operator with respect to the one in Medellín; (ii) in terms of energy losses, it is possible to observe a greater reduction in the urban area (11.76% more than in the rural area), which means that Medellín’s network operator saves 914.61 kWh when compared to the one in Capurganá; and (iii) by comparing the reductions in CO${}_2$ emissions, a behavior similar to that of the operating costs is observed: it is possible to achieve a reduction of 7.69% more than in the urban area, which translates into a reduction of 2547.17 kg of CO${}_2$ compared to the urban area.

These results are due to the fact that it is approximately 2.24 times more expensive to generate energy in Capurganá than it is in Medellín (Table 5), as the former is located in a non-interconnected electrical zone that requires a greater investment to transport diesel fuel for power generation. Similarly, when using diesel as fuel in rural areas, 1.62 times more kg of CO${}_2$ are emitted (Table 5). Therefore, it is possible to conclude that in rural areas, a greater economic and environmental impact is obtained in comparison with urban areas. However, due to the fact that in Capurganá, electricity is used for about 19 h a day on average [71], a greater technical impact is obtained in the city of Medellín, since this area uses electricity 24 h a day and at least 61% of the installed power is required (Figure 4).

One of the main barriers faced in the implementation of the proposed EMS is socio-economic in nature. In the country’s ZNIs, as is the case of Capurganá, local communities face difficulties due to: (i) their economic activities, such as fishing and agriculture; (ii) their low energy consumption and access to constant and high-quality energy; (iii) their high level of unsatisfied basic needs, such as the access to education; and (iv) their low payment capacity, as they do not have the necessary knowledge and capital to adopt, operate, and maintain PV generation systems, which is why the Colombian government promotes initiatives (such as this research paper) that seek to propose reliable energy alternatives and solutions based on the integration and management of non-conventional energy generation sources. Similarly, due to their complex cultural systems and their own worldviews, the availability of space for the installation of EMS is a complex issue, which can negatively affect the projects and investments made by the Colombian government. Therefore, it is important to educate the community of Capurganá about the operation of PV generators and the rational and efficient use of energy in order to improve the quality of life of the inhabitants of this ZNI by means of an EMS.

7. Conclusions and Future Work

The problem regarding the operation of PV generators in distribution grids was addressed in this paper by implementing an EMS based on a master–slave methodology. It is important to note the following. (i) The proposed EMS guarantees that the solutions obtained observe the technical-operating constraints of both test systems for the three simulation scenarios, as the nodal voltages are within the limits assigned for voltage regulation in each period of time. In the same way, the current that circulates through each conductor of the system in each period of time does not exceed the current thermal limit. (ii) The developed EMS allowed for an optimal dispatch of the solar resource by the PV generators, which can be implemented by any network operator, regardless of the geographical location of the electrical system, since the employed mathematical model contemplates the behavior curves of the solar resource and power consumption by users. Additionally, it was possible to show that the implementation of an EMS has greater economic and environmental impacts in distribution grids located in rural areas, and it has a greater technical impact in distribution grids located in urban areas.

On the other hand, numerical results for the urban (33-node test system) and rural (27-node test system) areas allow drawing the following conclusions:

✓

In one day of operation, for the urban area, a reduction of 2711.57 USD is obtained for the economic indicator, 1047.56 kWh for the technical indicator, and 3472.28 kg of CO${}_2$ for the environmental indicator, which represents reductions of 27.30, 31, and 27.69%, respectively. On the other hand, in rural areas, a reduction of 6521.44 USD is obtained for the economic indicator, 132.95 kWh for the technical indicator, and 6019.45 kg of CO${}_2$ for the environmental indicator, i.e., reductions of 35.16, 19.24, and 35.38%, respectively.

✓

In the urban area, a standard deviation of 0.0088% is obtained for the total operating costs, 0.0009% for the energy losses, and 0.0075% for the CO${}_2$ emissions, which is at least 200% lower in comparison with the other solution methodologies. As for the rural area, a standard deviation of 0.0007% is obtained for the total operating costs, 0.0005% for the energy losses, and 0.0004% for the CO${}_2$ emissions, i.e., at least 3000% lower with respect to the other solution methodologies.

✓

In the urban test system, the computation times are approximately 141.9966 s, 140.1495 s, and 140.5316 s for the three proposed simulation scenarios, respectively; while the rural test system reports about 125.9021 s, 126.52 s, and 131.3081 s. This demonstrates that the developed methodology allows solving multidimensional optimization problems with a continuous solution space (infinite combinations of power generation) at a low computational cost (less than 2.5 min) while guaranteeing the best response when compared to other metaheuristic algorithms.

As future work, it will be possible to solve the studied problem using metaheuristic algorithms with high numerical performance, e.g., the multi-verse optimizer, the salp swarm optimization algorithm, or the crow search algorithm. Likewise, it will be possible to consider the use of a multi-objective optimization approach that allows improving the economic, technical, and environmental indicators of the distribution grids in both rural and urban areas while observing their technical-operating conditions. An additional work could be extending the proposed formulation to distribution grid planning, including the problem regarding the integration of PV generators, as well as their investment costs and their carbon footprint throughout their useful life.