An Efﬁcient Methodology for Locating and Sizing PV Generators in Radial Distribution Networks Using a Mixed-Integer Conic Relaxation

: This paper proposes a new solution methodology based on a mixed-integer conic formulation to locate and size photovoltaic (PV) generation units in AC distribution networks with a radial structure. The objective function comprises the annual expected energy costs of the conventional substation in addition to the investment and operating costs of PV sources. The original optimization model that represents this problem belongs to the family of mixed-integer nonlinear programming (MINLP); however, the complexity of the power balance constraints make it difﬁcult to ﬁnd the global optimum. In order to improve the quality of the optimization model, a mixed-integer conic (MIC) formulation is proposed in this research in order to represent the studied problem. Numerical results in two test feeders composed of 33 and 69 nodes demonstrate the effectiveness of the proposed MIC model when compared to multiple metaheuristic optimizers such as the Chu and Beasley Genetic Algorithm, the Newton Metaheuristic Algorithm, the Vortex Search Algorithm, the Gradient-Based Metaheuristic Optimization Algorithm, and the Arithmetic Optimization Algorithm, among others. The ﬁnal results obtained with the MIC model show improvements greater than USD 100,000 per year of operation. All simulations were run in the MATLAB programming environment, using its own scripts for all the metaheuristic algorithms and the disciplined convex tool known as CVX with the Gurobi solver in order to solve the proposed MIC model.


General Context
Electrical distribution networks have undergone significant transformations due to the advancements made in power electronics, which allows for the optimal interface of renewable energy resources, energy storage systems, and controllable loads [1,2]. With the massive integration of these new actors into the electricity distribution sector, conventional passive grids have been transformed into active networks with multiple new challenges [3], such as minimizing the carbon footprint produced by the combustion of fossil fuels in thermal plants [4], improving grid efficiency by minimizing energy losses and by improving the voltage profiles [5], and reducing bad quality indicators regarding the reliability of distribution systems [6], among others.
One of the most important agents in the transformation of distribution networks corresponds to the important developments in small-scale renewable energy technologies that have changed the electricity consumption habits of millions of users in residential, industrial, and commercial applications at medium-and low-voltage levels [7][8][9]. These small generators can be massively integrated with distribution, and they can drastically change the behavior of the electrical network in terms of voltage profiles, energy losses, power flow directions, and the coordination of protective devices [10]. For this reason, several countries around the world have regulated the integration of small-scale generators in distribution networks to maintain economic and technical equilibrium in the distribution sector [11,12]. On the other hand, distribution companies can also integrate medium-size distributed energy resources in their distribution grids in order to minimize their carbon footprint and increase their profit from the distribution and commercialization of energy as a public service [13].
Regardless of the motivation for massively including renewable energy resources in electrical distribution networks, it is a fact that conventional passive distribution systems have stayed behind, and active networks require new research and developments in order to optimize their technical, economic, digital, and environmental dimensions [14].

Motivation
This research is motivated by the new challenges posed by the Sustainable Development Goals, one of them being the development of energy supply strategies for urban or rural users with zero greenhouse gas emitters [15]. To reach this objective, electric distribution companies, regulatory entities, and governments see renewable generation as a great opportunity to replace conventional thermal generation plants with renewable energy resources [16][17][18]. Solar photovoltaic (PV) and wind power systems are the most established, reliable, and mature technologies that can help with the continuous reduction in greenhouse gas emissions [19][20][21]; however, in countries located in the equatorial area (between the Capricorn and Cancer tropics), PV generation is preferred over wind power sources with regard to its integration into power systems at all voltage levels, given that the expected variations in solar resources throughout the year are minimal [22,23].

Literature Review
The problem regarding the optimal siting and sizing of renewable generation sources in electrical distribution networks has been largely explored through multiple optimization methodologies in the current literature. Some of the most recent reports are presented below.
The authors of [24] presented the application of a constructive heuristic algorithm based on the simulated annealing algorithm to locate wind power and solar PV sources as well as batteries in electrical distribution networks. Once the location of these devices was defined, a mixed-integer linear programming model was proposed to determine the optimal sizes for the renewables and batteries. Numerical results in three test feeders with 11, 135, and 230 nodes confirmed the effectiveness of the MINLP model at defining the size of the devices; however, the authors did not make any comparison with metaheuristic optimizers in order to confirm the effectiveness of the siting stage. In [25], the authors proposed a fuzzy-based methodology to determine the optimal regions for constructing PV generation units in order to deal with greenhouse gas emissions while considering uncertainties in the solar energetic resource. A validation of the proposed methodology in southern Iran demonstrated the effectiveness of combining the fuzzy-based approach with the analytical hierarchy process and Dempster-Shafer methods, in contrast with separately implementing these methodologies. Ref. [26] presented the application of the classical Chu and Beasley Genetic Algorithm (CBGA) with discrete-continuous codification in order to determine the optimal sizes of the PV generation units in electrical distribution networks. The objective of this research was to minimize the total annual grid operating costs while considering the energy purchasing costs in terminals of the substation and the investment and maintenance costs of the PV generation units installed. The solution of the exact MINLP model with the CBGA yielded better results for the IEEE 33-and IEEE 69-bus networks when compared to the implementation of the MINLP model in the General Algebraic Modeling System (GAMS) software and its BONMIN nonlinear solver. The authors of [15] discussed the main aspects of the massive integration of renewable power in power systems and their dependence on the simultaneous use of large-scale battery energy storage systems to make the power system reliable. A complete qualitative analysis regarding the intermittency of renewables and their percentage of penetration in the power system was provided with the help of several real scenarios around the world.
The problem regarding the optimal siting and sizing of PV generation units in distribution networks has been explored by [27]. The authors proposed the application of the Newton-Metaheuristic Algorithm (NMA) with a discrete-continuous codification. Numerical results in the IEEE 34-and IEEE 85-bus networks demonstrated the effectiveness of the NMA when compared to the CBGA and the BONMIN solver. A complete literature review regarding the optimal design of solar PV applications for rural electrification was presented by [28]. The authors of this research discussed the advantages of using smallscale solar systems to provide self-generation with the help of batteries as an alternative for improving the quality of life in regions without access to conventional power systems. Ref. [29] proposed a new combinatorial optimization method based on the modification of the gradient-based metaheuristic optimizer to locate and size PV sources in distribution networks. Its objective function is the minimization of the annual energy purchasing costs in terminals of the substation, as well as the installation and maintenance costs of the PV sources. The authors of [30] presented the comparison between photovoltaic and wind power systems in sub-transmission and distribution networks. Computational validations using a genetic algorithm in two test feeders composed of 9 and 33 buses show that both generation technologies efficiently improve the grid quality by means of their active and reactive power injection capabilities. Ref. [31] recently presented the application of the discrete-continuous vortex search algorithm (DCVSA) to locate and size PV sources in distribution networks with AC and DC topologies. The objective function under consideration involved the energy, investment, and installation costs of the PV generation units. Numerical results in the IEEE 33-and IEEE 69-bus networks demonstrated the superior performance of the DCVSA when compared to the BONMIN and the CBGA. Additional works that have addressed the problem regarding the optimal siting and sizing of PV sources in distribution networks include the Modified Arithmetic Optimization Algorithm (MAOA) [32], the Generalized Normal Distribution Optimization (GNDO) Algorithm [33], and the Tabu Search Algorithm [34], among others.
The main characteristic of the aforementioned optimization algorithms is that most of them belong to the family of combinatorial optimizers, which work with a master-slave (leader-follower) optimization design, where the master stage defines the optimal location and sizes of the PV sources and the slave stage is entrusted with the technical evaluation of this solution, i.e., the solution of the power flow problem.
It is worth mentioning that one of our proposal's distinct features is its solution approach, as this research focuses on the mathematical structure of the exact MINLP model that represents the studied problem, as well as the possibility of transforming it into a mixed-integer convex model [35,36].

Contribution And Scope
Considering the state-of-the-art presented above, this research makes the following contributions to the scientific literature: i.
The reformulation of the exact MINLP model to represent the problem regarding the optimal siting and sizing of PV sources in radial AC distribution networks through a Mixed-Integer Conic (MIC) model. The main advantage of the proposed MIC model is its solvability, which can ensure that the global optimum is reached via a combination of the Branch and Cut method with the interior point approach.
ii. Improved results in the IEEE 33-and IEEE 69-bus grids: about USD 100, 000/year per test feeder with respect to the best solution reported in the current literature with the GNDO algorithm [33].
It is important to mention that the results obtained in this research confirm that using maximum power point tracking for optimal dispatch is not the best option for operating PV generation units in electrical networks. With maximum power tracking, PV systems are forced to provide all their power to the grid at each period of time, whereas, with the possibility of optimal dispatch at each period of time, the distribution system can reach a better operational point. These observations are confirmed in the results section.
Note that the studied problem belongs to the family of distribution planning problems. In this study, we consider a planning period of 20 years for the evaluation of the NPV (net present value); however, the optimization model is solved for a typical daily operation scenario in the metropolitan area of Medellín, Colombia. In addition, it is important to mention that, once the PV generators are optimally located and sized with the proposed MIC model, it can be used as a conic dispatch model to determine changes in power consumption and renewable generation availability regarding the optimal daily operation of distributed energy resources in distribution networks by ensuring that the global optimum is reached [37,38].

Document Structure
The remainder of this research article is organized as follows: Section 2 describes the exact mathematical formulation that represents the problem regarding the optimal siting and sizing of PV generation units in radial AC distribution networks; Section 3 presents the mixed-integer conic reformulation based on the branch convex power flow formulation presented in [37]; Section 4 describes the main features of the IEEE 33-and IEEE 69-bus networks employed in the numerical validations, analysis, and discussions provided in Section 5; conclusions and proposals for future work are presented in Section 6.

Exact MINLP Model
The siting and sizing of PV generation units in electrical distribution networks with a radial structure implies a complex MINLP model given the combination of discrete and continuous variables [39], where i.
The integer and/or binary variables represent the nodes where the PV generation units will be placed; ii. The continuous variables are associated with the electrical variables such as the voltage and current magnitudes, the active and reactive power flows, and the active power generation in PV sources for each period of time, among others.
In addition, the nonlinear programming structure of the model is due to the nonlinear relation between voltages, currents, and powers in the active and reactive power balance constraints per node [37].
For the sake of simplicity, in the mathematical formulations presented below, lowercase letters correspond to variables, and uppercase letters represent constant values and parameters. In addition, it is worth mentioning that, in the proposed mathematical formulations, it is considered that the distribution networks have a strictly radial structure, and these are represented through a single-phase equivalent.

Objective Function Structure
To locate and size PV generators in radial AC distribution networks, an economical objective function is used as an objective function indicator, which has three main components: (i) the annual expected energy purchasing costs in terminals of the electrical substation, (ii) the investment costs regarding the installation of the PV generators, and (iii) the operating costs of the PV sources during the planning period. The formulation of this objective function and its components is presented below.
where z cost is the objective function value (net present value) of the annual energy purchasing cost in the substation (z a ), added to the investment costs in PV sources (z b ) and the operating and maintenance costs of the PV sources (z c ). C kWh is the expected average energy costs per year, T is the number of days in a year, R a is the expected return rate of the investments made by the distribution company, Y is the number of years of the planning period, R e is the expected rate of increment for the energy costs, and p sub si,h is the total power generation at the substation bus for each period of time h (power sent from node s to node i), where ∆h is defined as one hour (1 h). C pv is the cost of acquisition of a kW of PV power, p pv i is the nominal size of a PV generator connected at node i, C O&M is the maintenance and operation costs per kWh of energy generation in the PV sources, and p pv i,h is active power generation in the PV source connected at node i for the period of time h. The sets that contain all the nodes, the periods of time for daily operation, and the periods of time of the planning stage are denoted as N , H, and T , respectively. Remark 1. The main advantage of the objective function presented in Equation (1), with its components defined from (2) to (4), is that it is a convex objective function when the power generation at the substation bus (p sub si,h ) and the nominal sizes of the PV generation units are considered as variables, i.e., p pv i , taking into account that the remaining components are defined as constant parameters. Note that this objective function can be used with the same structure in the mixed-integer conic model presented in Section 3 due to its intrinsically convex structure.

Constraints
To locate and size PV generators in radial AC distribution networks, it is necessary to consider multiple linear and nonlinear constraints associated with the technical operation of the distribution network. These constraints include active and reactive power balance, voltage drops in lines, power transference per line, and son on. The complete list of constraints is presented below.
where p ij,h and q ij,h are the active and reactive power flow leaving node i towards destination j at time h; p jk,h and q jk,h are the active and reactive power sent from node j to any node k different from node i for each period of time; P d j,h and Q d j,h correspond to the total constant power consumption at node j for the period of time h; p pv j,h is the total power injected by the PV source connected at node j for each period of time; R ij and X ij represent the resistance and inductance effects in the distribution line in route ij; i ij,h is the variable associated with the current flow through the line that connects nodes i and j at each period of time; v i,h and v j,h correspond to the voltage magnitudes at nodes i and j for each period of time h, respectively; x j is a binary variable associated with the possibility of assigning a PV source at node j (x j = 1) or not (x j = 0); N max pv is a constant parameter related to the number of PV sources available for inclusion in the distribution network. The set that contains all the branches of the network is L.
Note that the set of constraints (5)-(13) can be understood as follows: Equations (5) and (6) correspond to the power balance equilibrium evaluated at each node per period of time; Equation (7) represents the voltage drop at each line, which is a function of the active and reactive power flows, as well as the branch parameters (i.e., resistance and inductance) and the current flow though the line that connects nodes i and j; Equation (8) corresponds to the application of Tellegen's second theorem to each sending node i in order to define the apparent power flow leaving it [40]; box-type constraint (9) defines the upper and lower active power generation bounds assigned for a PV source connected at node j for each period of time in the case that the binary variable x j is activated; box-type constraints (10) and (11) define the voltage regulation bounds applicable to all nodes and periods of time, as well as the thermal bounds regarding the distribution lines' current transport capacity; inequality constraint (12) defines the maximum number of renewable sources based on PV technology that can be included into the distribution network under study; Equation (13) describes the binary nature of the decision variable x j . Figure 1 summarizes the main characteristics of the optimization model (1)- (13). Note that the studied optimization model is only nonlinear and non-convex in its continuous part, which is associated with the active and reactive power balance constraints, voltage drops per line, and the calculation of the apparent power [13].

Remark 2.
As demonstrated in [37], it is possible to transform the set of nonlinear constraints (5)-(8) into a conic set of constraints, which allows transforming the exact MINLP model (1)-(13) into a mixed-integer conic (MIC) model, as presented in the next section.

Mic Reformulation
The main contribution of this research is the reformulation of the exact MINLP model (1)- (13), especially the subset of constraints (5)-(8) as a set of conic constraints based on the branch flow theory for distribution networks presented in [37]. To this effect, the following set of auxiliary variables is used: [41]. Note that, with the aforementioned definition of the auxiliary variables, the set of constraints (5)-(7) can be directly transformed into a set of linear affine (i.e., convex) constraints. These equations take the following form: To become a conic constraint, Equation (8) uses the concept of hyperbolic equivalence between the product of two variables [37].
Equation system (17) is still non-linear and non-convex due to the use of the equality symbol since it corresponds to the surface of a hypersphere, which requires additional relaxation to become a conic convex constraint [37]. To achieve this, the equality symbol is replaced with a lower-equal symbol, as defined in Inequality (18).
A modification of the objective function (1) is made in order to consider the effect of the annual costs of the energy losses while deciding the optimal siting and sizing of the PV sources. This effect is only included in the conic model to improve its solution properties, since power losses can be represented as quadratic (as a function of the variable i ij,h ) or as linear (as a function of the auxiliary variable l ij,h ) function: Note that β is the factor assigned to the calculation of the energy losses. In order to control its effect on the objective function value, this factor takes three possible values: 0, 1/2, and 1. The location and size of the PV source with each factor is illustrated in the results section. Figure 2 illustrates the transformation of the MINLP model presented in Figure 1, where there are only linear affine equations and binary constraints.

Remark 3.
To ensure that the global optimum is found via the proposed MIC model, it is possible to implement the Branch and Cut method in conjunction with the interior point approach, as presented in [42], to deal with mixed-integer conic models.
In order to guarantee that the solution is reached by the MIC model, it is evaluated in the exact MINLP model by fixing all the binary variables on it, which reduces it to a daily optimal power flow model, thus ensuring that the exact behavior of the objective function is known, as well as the PV generation outputs, as recommended in [13].

Test Feeders and Model Parameters
This section presents the main characteristics of the test feeders employed for all the numerical validations, as well as a model characterization to calculate the objective function and the behavior of the demand and solar generation curves.

IEEE 33-Bus Network
This network is a radial electrical distribution system composed of 33 nodes and 32 lines, which is operated with a line-to-ground voltage magnitude of 12.66 kV in terminals of the substation connected at bus 1. The electrical configuration of this test feeder is depicted in Figure 3a, and its electrical data are reported in Table 1.   Node i Node j R ij (Ω) X ij (Ω) P j (kW) Q j (kvar) Node i Node j R ij (Ω) X ij (Ω) P j (kW) Q j (kvar)

IEEE 69-Bus Network
This system is also a radial distribution network with 69 nodes and 68 lines operated with a line-to-ground voltage of 12.66 kV in terminals of the substation connected at node 1. The electrical configuration of this test feeder is depicted in Figure 3b, and its electrical parameters are presented in Table 2.

Parameters for the Economic Assessment
To evaluate the objective function and its components, i.e., Equations (1)-(4), all the parameters listed in Table 3 are considered. The mean average generation and demand curves for the metropolitan area of Medellín, Colombia, were used to model the daily behavior of the demand and generators in two test feeders. These curves are plotted in

Computational Validations
The proposed MIC model for siting and size PV generation units in AC distribution networks with a radial structure was implemented in the MATLAB software version 2021b from MathWorks (Natick, MA, USA). To solve the MIC model, the convex disciplined tool known as CVX was used with the Gurobi solver. All the implementations were conducted on a desktop computer with an Intel(R) Core(TM) i7-7700 2.8-GHz processor and 16.0 GB of RAM running a 64-bit version of Microsoft Windows 10 Home.

Remark 4.
The solutions reached with each metaheuristic optimization algorithm assumed that the PV generation units were operated using the maximum power point tracking algorithm, i.e., the PV curve in Figure 4 was perfectly followed; however, in the case of the proposed MIC model, the PV generation units were freely operated, i.e., they did not necessarily reach the maximum power depicted in the PV curve of Figure 4.
In order to implement all the comparative metaheuristic optimizers, a discrete-continuous codification based on the proposal presented in [26] is employed. This codification takes the structure defined in (20).
where X sol represents a possible solution individual in which the first N max pv components are related to discrete variables that define the potential nodes where the PV sources will be assigned, and the positions N max pv + 1 to 2N max pv determine the potential sizes assigned to these PV sources.
Due to the fact that metaheuristic optimizers such as the gradient-based or the Newtonbased approaches are defined in the continuous domain, in order to apply their evolution rules, the integer nature of the decision variables is relaxed; however, once the new individual is found, it is rounded in its first N max pv to maintain the feasibility of the solution space during the exploration and exploitation searching processes. Table 4 presents the numerical performance of all the metaheuristic optimizers and the proposed MIC model for the IEEE 33-bus network. The numerical results in Table 4 show that:

IEEE 33-Bus Network
i. The best metaheuristic approach for the IEEE 33-bus grid was the GNDO method, with an objective function of USD/year 2, 699, 671.75; and the worst metaheuristic optimizer was the NMA, with an objective function value of USD/year 2, 700, 227.33; however, both solutions only differed by about USD/year 555.58. These values imply that all solution methodologies based on the application of metaheuristic methods reported in Table 4 are contained between both bounds. ii. The proposed MIC model with different β-values found objective function values between USD/year 2, 603, 465.00 and USD/year 2, 597, 139.00, i.e., there is a difference of about USD/year 6326 when β = 0 and β = 1. iii. By comparing the GNDO solution (the best solution among the metaheuristic approaches) and the MIC solution with β = 1, the difference was about USD/year 102, 532.75. This confirms that the operation of the PV generators without maximum power point tracking is considerably better than the approaches where the PV sources are forced to follow the maximum power point (see metaheuristic results). iv. As for the places and sizes of the PV generators, it was observed that, in the case of the MIC model, the β−parameter (i.e., the effect of the power losses) has a significant impact on the final nodal location of the PV generation sources, and well as on their sizes; however, in the objective function, this effect is minimized due to the multimodal behavior of the solution space in the exact MINLP formulation. This behavior of the studied problem also explains the multiplicity in the set of solutions reported by all the metaheuristic optimizers.
It is worth mentioning that, regarding processing times, the proposed MIC models for the IEEE 33-bus grids take about 45 s to solve the conic model, which is an adequate processing time considering the dimension of the solution space, as there are about 4960 options for discrete variables, with the main advantage that the final solution reached is indeed the global optimal one.
Regarding the annual expected improvements in the objective function, Figure 5 presents the percentage of reductions reached by each optimization method with respect to the benchmark case.
The behavior of the expected annual reductions implies the following: (i) all metaheuristic optimizers that operated using the concept of maximum power point tracking found solutions with reductions between 27.03 and 27.04%, i.e., a difference lower than 0.01%, which allows concluding that, in the IEEE 33-bus network, all the metaheuristic methods have a very good performance (however, these are indeed local optimal solutions, since it is not possible to prove that these algorithms reach the global optimal solution due to their random-based nature); (ii) the MIC model reached reductions between 29.64 and 29.82%, where the difference is associated with the inclusion of the β−coefficient in the objective function. Nevertheless, all of the MIC solutions improved the objective function values by more than USD 100, 000 with respect to the metaheuristic approaches. As mentioned above, this improvement is due to the fact that the MIC model allows operating the PV sources not necessarily with maximum power point tracking. Note that each one of the MIC models provides an optimal solution, but they differ from each other because each one depends on the value assigned to the β parameter, which implies that each objective function is different.   Table 5 presents the comparative analysis between the performance of the metaheuristic optimizers reported in the specialized literature and that of the proposed MIC model with different β−coefficients. The numerical results in Table 5 show that: As for processing times, in the IEEE 69-bus grid, the proposed MIC models take about 438 s to solve the conic model, which is an adequate processing time considering the dimension of the solution space, as there are about 50116 options regarding discrete variables. Note that, even though this processing time is high when compared to metaheuristic approaches, it has the advantage that the global optimal solution is reached with 100% of certainty, which is not possible with any metaheuristic model due to their random-based nature.

IEEE 69-Bus Network
On the other hand, Figure 6 depicts the expected reductions of the metaheuristic and the MIC models in the IEEE 69-bus network when compared to the benchmark case.  In the case of the IEEE 69-bus network, all the analyzed metaheuristic methods reached improvements between 27.12 and 27.16% regarding the objective function, i.e., there were differences lower than 0.04%. This small difference confirmed that all the metaheuristic optimizers are adequate to solve the problem concerning the optimal placement and sizing of PV generation units; however, these solutions remain classified as sub-optimal due to the impossibility of theoretically ensuring their global nature. The MIC models showed additional improvements of about 1.88 and 2.67% in the objective function value when compared to the GNDO approach (the best metaheuristic method). These improvements confirmed that, when the maximum power point is not tracked, it is possible to reach better results regarding the final objective function value when studying the problem of the optimal placement and sizing of PV generators in distribution networks, with the main advantage that the solutions provided by the models are the global optimum for each assigned value in the β parameter.

Conclusions
The problem regarding the optimal siting and sizing of PV generation units in radial AC distribution networks was addressed in this research by transforming the exact MINLP model into a MIC one. Numerical results in the IEEE 33-and IEEE 69-bus networks showed that: (i) in both test feeders, the best metaheuristic optimization method was the GNDO method, with improvements in the final objective function value of about 27.04 and 27.16%, respectively; (ii) all the metaheuristic optimization methods showed improvements higher than 27.03% in the IEEE 33-bus network and 27.12% in the IEEE 69-bus network, which confirmed that all the studied combinatorial optimization approaches are suitable for solving the problem under study with excellent numerical results; (iii) operating the PV generators while applying maximum power point tracking showed solutions higher than USD/year 2, 699, 600.00 for the IEEE 33-bus system, while the IEEE 69-bus network showed solutions higher than USD/year 2, 824, 900.00; however, when the grids are freely operated during the day, these values can reduced by about USD/year 100, 000 per test feeder, as was demonstrated with the MIC model for different β−coefficients.
The main advantage of the proposed MIC model is that it is possible to ensure that the global optimum is found via the interior point method combined with the Branch and Cut optimization approach, which implies that the solution for each β−coefficient will be the same in each run; however, when metaheuristic optimization methods are applied to solve MINLP models, statistical tests are required in order to guarantee their adequate performance. Nevertheless, global optimum finding properties can not be attributed to these methodologies.
As future work, it will be possible to conduct the following studies: (i) combining PV generation units with reactive power compensators for annual cost reductions in distribution networks and (ii) extending the proposed MIC model to direct current networks with monopolar and bipolar configurations. Funding: This research was supported by Minciencias, Universidad Nacional de Colombia, Universidad del Valle, and Instituto Tecnológico Metropolitano under the research project "Dimensionamiento, planeación y control de sistemas eléctricos basados en fuentes renovables no convencionales, sistemas de almacenamiento y pilas de combustible para incrementar el acceso y la seguridad energética de poblaciones colombianas" (Minciencias code 70386), which belongs to the research program "Estrategias para el desarrollo de sistemas energéticos sostenibles, confiables, eficientes y accesibles para el futuro de Colombia" (Minciencias code 1150-852-70378, Hermes code 46771).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used in this study are reported in the paper's figures and tables. p jk,h Active power flow leaving node j towards destination k at time h (W). q jk,h Reactive power flow leaving node j towards destination k at time h (var). P d j,h Total active power consumption at node j for period of time h (W).

Q d j,h
Total reactive power consumption at node j for period of time h (var). p pv j,h Total power injected by the PV source connected at node j for each period of time h (W).

R ij
Resistance effect in the distribution line in route ij (Ω).

X ij
Reactance effect in the distribution line in route ij (Ω). i ij,h Current flow through the line that connects nodes i and j in each period of time h (A). v i,h Voltage magnitude at node i for each period of time h (V). v j,h Voltage magnitude at node j for each period of time h (V).
x j Binary variable associated with the possibility of assigning a PV source at node j (x j = 1) or not (x j = 0). β Factor that allows including or not the expected annual costs of the energy loss.