Efficient Integration of Photovoltaic Solar Generators in Monopolar DC Networks through a Convex Mixed-Integer Optimization Model

Abstract

The problem regarding the optimal siting and sizing of photovoltaic (PV) generation units in electrical distribution networks with monopolar direct current (DC) operation technology was addressed in this research by proposing a two-stage convex optimization (TSCO) approach. In the first stage, the exact mixed-integer nonlinear programming (MINLP) formulation was relaxed via mixed-integer linear programming, defining the nodes where the PV generation units must be placed. In the second stage, the optimal power flow problem associated with PV sizing was solved by approximating the exact nonlinear component of the MINLP model into a second-order cone programming equivalent. The main contribution of this research is the use of two approximations to efficiently solve the studied problem, by taking advantage of convex optimization models. The numerical results in the monopolar DC version of the IEEE 33-bus grid demonstrate the effectiveness of the proposed approach when compared to multiple combinatorial optimization methods. Two evaluations were conducted, to confirm the efficiency of the proposed optimization model. The first evaluation considered the IEEE 33-bus grid without current limitations in all distribution branches, to later compare it to different metaheuristic approaches (discrete versions of the Chu and Beasley genetic algorithm, the vortex search algorithm, and the generalized normal distribution optimizer); the second simulation included the thermal current limits in the model’s optimization. The numerical results showed that when the maximum point power tracking was not regarded as a decision-making criterion, the expected annual investment and operating costs exhibited better performances, i.e., additional reductions of about USD 100,000 in the simulation cases compared to the scenarios involving maximum power point tracking.

1. Introduction

Electrical energy generation is a costly process (both economically and environmentally) that increasingly consumes our natural resources [1,2,3]. In line with sustainable development goals, worldwide efforts are being made towards an energy transition from fossil sources to renewable and clean energy markets [4,5,6]. In the electrical sector, one of the main aspects of the energy transition is the use of medium- and large-scale renewable generation systems (mainly photovoltaic and wind sources) combined with energy storage systems, to replace—to the greatest possible extent—the participation of thermal sources (fossil-based sources) in the energy supply chain for all users of electrical energy [7,8,9,10].

Electrical distribution is one of the sectors of the electricity supply chain that is also undergoing constant change, as multiple governments and regulatory policies have driven the massive integration of renewable energy resources in these systems [11]. The main idea is to transform conventional passive distribution networks with unidirectional power flows into active distribution grids with a high participation of renewable energy resources, aiming to reduce the carbon footprint of conventional power systems at the medium- and low-voltage levels [12].

This research aims to contribute to the energy transition by proposing an optimization methodology to determine the best possible location and size of photovoltaic (PV) generation plants in direct current (DC) distribution networks with monopolar configurations. In the specialized literature, the optimal integration/operation of renewable energy plants in medium- and low-voltage grids has been widely studied. Some of the recent advances in this topic are presented below.

The authors of [13] presented an efficient technique for operating PV systems while considering the maximum power tracking point, in order to improve the efficiency, stability, and reliability of these generators. An adaptive neural fuzzy interference system was proposed, to design the maximum power tracking, given its fast tracking speed at reaching steady state and low oscillation. A standalone network with a large historical dataset was used to validate the proposed methodology via MATLAB/Simulink.

The work by [11] presented a pre-feasibility analysis regarding the integration of small-scale renewable energy resources based on wind and PV power under Colombian Law 1715. The authors explored the residential benefit of using renewable energy resources for self-energy supply as a function of the economic conditions of each residential user (using economic strata to classify all end users). Numerical validations demonstrated that investments could not be amortized during the lifetime of small wind turbines in any strata; however, investments in small PV plants could be amortized in between 6.2 and 8.7 years. In addition, for strata 4 to 6 (i.e., middle- to high-class income), these investments could be recovered in the fourth year of operation.

The authors of [14] applied the Jaya optimization algorithm, to locate and size PV plants in medium-voltage distribution networks. The numerical results in the IEEE 33-bus grid, regarding nodal location and PV sizing, yielded better results with regard to minimizing power losses and improving voltage profiles against well-known optimizers, such as genetic algorithms and particle swarm optimization, among others.

In [15], a two-stage optimization approach was proposed, to locate renewable generation sources and energy storage systems in medium- and low-voltage distribution networks. In the first stage, the simulated annealing optimization method defined the nodes where the batteries and renewable energy resources needed to be located. In the slave stage, a convex optimization strategy was proposed, to determine the size of the renewables and batteries, in order to minimize investment and operating costs during the planning period. Numerical results in different test feeders, with between 11 and 230 nodes, demonstrated the effectiveness of the linearized convex model when compared to second-order cone formulations.

The authors of [16] applied the generalized normal distribution optimizer, to determine the optimal location and sizing of PV generation sources in monopolar DC distribution networks via a master-slave optimization strategy. In the slave stage, the successive approximations power flow method was implemented, to evaluate the multi-period optimal power flow problem, aiming to define the annual operating costs of the grid, i.e., the energy production costs in the slack node, and the maintenance and operating costs of the PV plants. The master stage determined the PV location and sizes, using a discrete–continuous codification. The numerical results in the IEEE 33- and IEEE 69-bus grids demonstrated the effectiveness of the proposed optimization approach.

Additional literature approaches regarding the optimal integration (location/sizing and operation) of PV plants in electrical distribution networks include the following optimization methods: particle swarm optimization [17,18]; the water cycle algorithm [19]; the gravitational search algorithm [20]; the vortex search algorithm [21]; the crow search algorithm [22], and a mixed-integer programming approximation [23], among others.

Considering the literature, this research contributes to the following aspects: (i) a new two-stage optimization methodology that defines the nodes where the PV plants must be located, via a mixed-integer linear programming formulation based on simplifying the branch power flow model; in this first stage, the binary component of the exact mixed-integer linear programming (MINLP) model that defines the optimal PV plant placement and sizing problem is solved; (ii) once the binary part of the MINLP model has been solved, the resulting nonlinear programming model is transformed into its second-order cone programming equivalent, to define the optimal size of the PV plants at the nodes where they were previously located. This stage considers the possibility of operating the PV plants with maximum power point tracking or free hourly dispatch, with a high impact on the final numerical results.

The numerical results in the monopolar version of the IEEE 33-bus grid confirmed the effectiveness of the proposed two-stage convex optimization approach, in comparison to the solutions reported by the generalized normal distribution system optimizer, the Chu and Beasley genetic algorithm, and the vortex search algorithm [16], respectively. In this research, no uncertainties were considered in the PV generation and demand profiles, i.e., these data were assumed to be constant input parameters. However, for future research, the stochastic behavior of these variables could be considered [21]. In addition, the proposed two-stage convex optimization approach only applies to purely radial distribution networks. This is because the second-order cone programming model, based on the branch power flow formulation, was developed for radial distribution networks [24]. However, it can be considered an opportunity for research to propose an alternative approximated formulation approach that allows viewing radial and meshed grids, to maintain the convexity of the solution space. An important area of interest for future works would be to consider the demand and generation of variable curves throughout the planning period, to better approximate the net profits expected by the utility company. In the specialized literature, some authors have proposed applying stochastic optimization methods, to deal with generation and demand uncertainties, and these can be consulted in [25,26].

The remainder of this document is organized as follows. Section 2 presents the exact MINLP formulation of the problem regarding the optimal placement and sizing of PV sources in monopolar DC distribution networks. Section 3 describes the proposed two-stage optimization methodology, which is based on a mixed-integer linear programming approximation to determine the nodes where the PV sources must be located in the first optimization stage, while the second stage reformulates the nonlinear programming model, regarding the hourly optimal power flow, as a second-order conic programming model. Section 4 shows the main characteristics of the test feeder under analysis, which corresponds to the monopolar DC version of the IEEE 33-bus grid. Section 5 presents the main numerical results of the proposed two-stage optimization method, in comparison to combinatorial optimization methods, and their complete analysis and discussion. Finally, Section 6 lists the main conclusions of this research, and some possible future works.

2. Mathematical Formulation

The problem regarding the optimal design of PV systems in electrical networks is a challenging task in the areas of electrical analysis and optimization, given that:

i.

The electrical network itself is modeled as a set of nonlinear non-convex equations associated with Kirchhoff’s laws applied to electrical circuits with constant power terminals (i.e., nonlinear loads);

ii.

The integration of a group of shunt devices (renewable generation or batteries) introduces binary decision variables that transform the solution space into that of a disjoint optimization problem, i.e., a problem with a mixed-integer solution space.

The complete optimization model that represents the optimal siting (binary) and sizing (nonlinear and continuous) problems for PV generation units in monopolar DC distribution networks is detailed below.

2.1. Objective Function

The main objective of integrating dispersed generation in electrical networks is to minimize the annual grid operating costs: these costs can comprise multiple components, such as energy purchasing, the maintenance of PV plants, and energy losses. This research employed an objective function associated with energy purchasing and maintenance. The mathematical structure of the objective function is defined in Equations (1)–(3) [15,27].

$$\begin{array}{c}{A}_{\mathrm{cost}}=D_1+D_2,\end{array}$$

(1)

$$\begin{array}{c}{D}_1=C_{\mathrm{kWh}}T\left(\frac{t_a}{1-{\left(1+t_a\right)}^{-N_t}}\right)\left(\sum _{h\in \mathcal{H}}\sum _{i\in \mathcal{N}}{P}_{0i,h}^{cg}\Delta h\right)\left(\sum _{t\in \mathcal{T}}{\left(\frac{1+t_e}{1+t_a}\right)}^t\right)\end{array}$$

(2)

$$\begin{array}{c}{D}_2=C_{\mathrm{pv}}\left(\frac{t_a}{1-{\left(1+t_a\right)}^{-N_t}}\right)\left(\sum _{i\in \mathcal{N}}{P}_i^{pv}\right)+C_{\mathrm{O}\&\mathrm{M}}T\left(\sum _{h\in \mathcal{H}}\sum _{i\in \mathcal{N}}{P}_{i,h}^{pv}\Delta h\right).\end{array}$$

(3)

In Equation (1), the component $A_{\mathrm{cost}}$ describes the expected annual operating costs of the network, which are the sum of the components $D_1$ and $D_2$ (see Equations (2) and (3)), with the former term being the expected energy purchasing costs at the terminals of the substation, and the latter being the expected maintenance and operating costs of the PV generation units. In addition, $C_{\mathrm{kWh}}$ corresponds to the average value of a kilowatt of power generated in the substation bus; T is a constant associated with the number of days in a year; $t_{\mathrm{a}}$ is the expected return rate of the investments made by the distribution company; $N_t$ is the number of years in the planning period; $P_{0i,h}^{cg}$ corresponds to the power generation at the slack node (substation bus), which is defined as the power flowing from the substation (node 0) to any node i interconnected with it; $\Delta h$ represents the time variation regarding the daily operation intervals (1 h); $t_{\mathrm{e}}$ means the expected rate of increasing the energy generation costs per year of operation; $C_{\mathrm{pv}}$ denotes a parameter associated with the investment costs in a kWp of PV-based power; $P_i^{pv}$ is the expected size of a PV generation installed at node i, which generates $P_{i,h}^{pv}$ per period of analysis; and $C_{\mathrm{O}\&\mathrm{M}}$ represents the maintenance and operating costs coefficient per kWh of energy generated by a PV generation unit. In addition, $\mathcal{H}$ is the set that contains all the periods in the daily operating scenario, and $\mathcal{N}$ is the set associated with the number of nodes of the monopolar DC network under analysis.

Remark 1. 

Note that the mathematical structure of the objective function in (1) is a linear combination of components (2) and (3), which have a convex structure, as both are linear functions of the decision variables $P_{0i,h}^{cg}$, $P_i^{pv}$, and $P_i^{pv}$. This objective function property was used to propose the two-stage optimization approach to locating and sizing PV generation units in monopolar DC networks.

2.2. Set of Constraints

Multiple technical constraints regarding Kirchhoff’s laws must be satisfied, in order to ensure the optimal operation of the monopolar DC network after inserting the PV generation units.

The power equilibrium constraint corresponds to applying Kirchhoff’s first law while using a power structure. This equality constraint is defined in (4) [28]. Note that, in order to obtain the power balance constraint in (4), the branch power flow formulation for direct current networks is employed; further details can be consulted in [29].

$$\begin{array}{c}{P}_{ij,h}-R_{ij}{i}_{ij,h}^2\sum _{k:\left(jk\right)\in \mathcal{L}}{P}_{jk,h}=P_{j,h}^d-P_{j,h}^{pv},\left\{\begin{array}{c}\forall j\in \mathcal{N},j\ne slack,\\ \forall h\in \mathcal{H}\end{array}\right\},\end{array}$$

(4)

where $P_{ij,h}$ ($P_{jk,h}$) denotes the power flowing from node i (j) to node j (k) per period h; $R_{ij}$ corresponds to the resistive parameter associated with the conductor assigned to the route $ij$; $i_{ij,h}$ represents the current flow in the route $ij$ at each period; and $P_{j,h}^d$ is the power demanded at node j in the period h. Note that $\mathcal{L}$ is the set that defines the routes (distribution lines) in the studied test feeder.

To associate the power balance constraint in (4) with the variation in the voltage variables, the voltage drop in each distribution line is represented as a function of its resistance and the power flow, as defined in (5) [24]:

$$\begin{array}{c}{v}_{j,h}^2=v_{i,h}^2-2R_{ij}{P}_{ij,h}+R_{ij}^2{i}_{ij,h}^2,\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\},\end{array}$$

(5)

where $v_{j,h}$ ($v_{i,h}$) is the voltage variable at node j (i) per period of analysis.

To relate voltage variables per node to the current variables, Tellegen’s second theorem is applied, as presented in (6) [29]:

$$\begin{array}{c}{P}_{ij,h}=v_{i,h}{i}_{ij,h},\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\}.\end{array}$$

(6)

To determine the expected size of the PV generation plant to be connected at node i, the following binary constraint is implemented [30]:

$$\begin{array}{c}{z}_j{P}_{pv}^{min}\le P_j^{pv}\le z_j{P}_{pv}^{max},\left\{\forall j\in \mathcal{N}\right\},\end{array}$$

(7)

where $z_j$ is the binary variable that allows for defining whether the PV generator is located at node j ($z_j=1$) or not ($z_j=0$), and $P_{pv}^{min}$ and $P_{pv}^{max}$ represent the minimum and maximum sizes allowed for the PV generation units.

The following inequality constraint is added to the optimization model, in order to determine the hourly expected power generation by each PV source integrated into the monopolar DC network [21]:

$$\begin{array}{c}0\le P_{j,h}^{pv}\le G_h^{pv}{P}_j^{pv},\left\{\forall j\in \mathcal{N},\forall h\in \mathcal{H}\right\},\end{array}$$

(8)

where $G_h^{pv}$ denotes the expected profile of the energy generation curve for the PV generation plants, which is provided as a prediction by the distribution company [31,32].

To ensure that the voltage profiles in all the network nodes fulfill the regulatory policies, the optimization model includes the box-type constraint (9). According to this constraint, the voltage magnitudes must be between their minimum and maximum values (i.e., $v^{min}$ and $v^{max}$) [15]:

$$\begin{array}{c}{v}^{min}\le v_{j,h}\le v^{max},\left\{\forall j\in \mathcal{N},\forall h\in \mathcal{H}\right\}.\end{array}$$

(9)

To protect the distribution branches, with regard to their thermal current limitations, box-type constraint (10) ensures that the current flowing through each distribution line is within its limits [30]:

$$\begin{array}{c}-i_{ij}^{max}\le i_{ij,h}\le i_{ij}^{max},\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\},\end{array}$$

(10)

where $i_{ij}^{max}$ is the thermal bound associated with the conductor connected in the $i{j}^{th}$ route.

Finally, the maximum number of PV generators available for installation (i.e., $N_{pv}^{max}$) is set as an inequality constraint in (11). In addition, the binary nature of the decision variable is ratified in (12) [21]:

$$\begin{array}{c}\sum _{j\in \mathcal{N}}{z}_j\le N_{pv}^{max},\end{array}$$

(11)

$$\begin{array}{c}{z}_j\in (0,1),\left\{\forall i\in \mathcal{N}\right\}.\end{array}$$

(12)

Remark 2. 

Note that, in Set (4)–(12), most constraints (66.67%) define a mixed-integer convex solution space, as is the case in (8)–(12); however, three non-convex constraints must be analyzed and redefined, to reach a complete convex optimization model, i.e., (4)–(6).

Due to constraints (4)–(6), the complete optimization model (1)–(12) belongs to the family of MINLP formulations, which implies that there are still challenges regarding its optimal solution: this constitutes the contribution of this research, as it proposes a new approximated solution strategy based on two convex approximations, in order to provide an efficient solution to the exact MINLP model with regard to the siting and sizing of PV generation units in monopolar DC networks. In addition, it is essential to highlight that, as for the solution space, this research assumed that it is possible to install a set of PV generation units in each node without considering any physical constraint regarding the physical area available for these installations: this is a typical assumption for locating dispersed generation units in electrical networks, given that it is the worst possible scenario, as it exhibits the largest possible dimensions of the solution space [33].

3. Proposed Two-Stage Optimization Methodology

This section presents the proposed two-stage optimization approach to locating and sizing PV generators in monopolar DC networks, with the purpose of minimizing the expected annual operating costs. The first stage corresponds to a relaxation of the optimization model, by neglecting the effect of the voltage profile, in order to obtain a mixed-integer linear programming model that allows for identifying the set of nodes where the PV sources must be located. The second stage uses a mixed-integer conic approximation, to determine the expected size of the PV generators, while considering the binary variables to be fixed.

3.1. Node Selection for the PV Generation Units

To obtain a mixed-integer convex approach that allowed for selecting the set of nodes where the PV generators needed to be installed, this study used an auxiliary variable associated with the current square value, i.e., $l_{ij,h}=i_{ij,h}^2$. With this auxiliary variable, the power balance constraint (4) took the following form:

$$\begin{array}{c}{P}_{ij,h}-R_{ij}{l}_{ij,h}\sum _{k:\left(jk\right)\in \mathcal{L}}{P}_{jk,h}=P_{j,h}^d-P_{j,h}^{pv},\left\{\begin{array}{c}\forall j\in \mathcal{N}\\ j\ne slack\\ \forall h\in \mathcal{H}\end{array}\right\}.\end{array}$$

(13)

Remark 3. 

To obtain a simplified convex model that allowed for determining the set of nodes where the PV generation units needed be located, it was assumed that the voltage variations were minimal with respect to the reference value, which implied that $v_{i,h}\approx v_{j,h}\approx 1$. With this assumption, the contribution of the equality constraint (5) to the optimization model could be neglected.

In light of the above-presented considerations on the current variable, a new term was added to the objective function, with the aim of emulating the effect of energy losses during the operation of the PV plants of a monopolar DC network. Note that, according to (6), $P_{ij,h}^2\approx i_{ij,h}^2=l_{ij,h}$, which implied that a component considering the effect of the energy losses could be added to the costs of said losses, as follows:

$$\begin{array}{c}\hfill min{z}_{\mathrm{approx}}=A_{\mathrm{costs}}+C_{kWh}T\sum _{h\in \mathcal{H}}\sum _{ij\in \mathcal{L}}{R}_{ij}{l}_{ij,h}{\Delta }_h.\end{array}$$

(14)

Note that the relaxed optimization model used to identify the set of nodes where the PV generation sources needed be installed took the form presented in (15).

$$\begin{array}{cc}\hfill \mathrm{Obj}.\mathrm{Func}.:& \hfill \\ \hfill & min{z}_{\mathrm{approx}}=A_{\mathrm{costs}}+C_{kWh}T\sum _{h\in \mathcal{H}}\sum _{ij\in \mathcal{L}}{R}_{ij}{l}_{ij,h}{\Delta }_h.\hfill \\ \hfill \mathrm{Subject}\mathrm{to}:& \hfill \\ \hfill & P_{ij,h}-R_{ij}{l}_{ij,h}\sum _{k:\left(jk\right)\in \mathcal{L}}{P}_{jk,h}=P_{j,h}^d-P_{j,h}^{pv},\left\{\begin{array}{c}\forall j\in \mathcal{N},\\ j\ne slack,\\ \forall h\in \mathcal{H}\end{array}\right\},\hfill \\ \hfill & P_{ij,h}=l_{ij,h},\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\}\hfill \\ \hfill & z_j{P}_{pv}^{min}\le P_j^{pv}\le z_j{P}_{pv}^{max},\left\{\forall j\in \mathcal{N}\right\},\hfill \\ \hfill & 0\le P_{j,h}^{pv}\le G_h^{pv}{P}_j^{pv},\left\{\forall j\in \mathcal{N},\forall h\in \mathcal{H}\right\},\hfill \\ \hfill & {\left(v^{min}\right)}^2\le u_{j,h}\le {\left(v^{max}\right)}^2,\left\{\forall j\in \mathcal{N},\forall h\in \mathcal{H}\right\},\hfill \\ \hfill & -{\left(i_{ij}^{min}\right)}^2\le l_{ij,h}\le {\left(i_{ij}^{max}\right)}^2,\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\},\hfill \\ \hfill & \sum _{j\in \mathcal{N}}{z}_j\le N_{pv}^{max},\hfill \\ \hfill & z_j\in (0,1),\left\{\forall i\in \mathcal{N}\right\}.\hfill \end{array}$$

(15)

The main characteristic of the optimization model (15) is that it belongs to the family of mixed-integer linear programming, which implies that with the branch and cut approach it is possible to obtain its optimal solution [34]. Finally, the result of interest in this model was the set of nodes where the PV plants must be placed (i.e., $z_j$), which was set as an input vector for the second stage, with regard to the optimization model that would define the set of PV plants.

3.2. Calculating the Size of the PV Generation Units

Once the mixed-integer linear programming model (15) had been solved, the location of the PV generators was contained in the variable $z_j$, which was considered as an input parameter in this optimization stage. To determine the expected size of the PV plants, the optimization model (1)–(12) was rewritten as a second-order cone programming model [24]. To obtain the conic model, this research introduced an additional auxiliary variable regarding voltages (i.e., $u_{j,h}=v_{j,h}^2$), which allowed redefining for constraint (5) as (16):

$$\begin{array}{c}{u}_{j,h}=u_{i,h}-2R_{ij}{P}_{ij,h}+R_{ij}^2{i}_{ij,h},\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\},\end{array}$$

(16)

which was an affine constraint, in terms of the new auxiliary variables $u_{i,h}$, $u_{j,h}$, and $l_{ij,h}$.

Now, the only non-convex constraint in the optimization model (1)–(12) corresponded to the definition of the power flow for each branch, i.e., Equation (6); however, this equation could be convexified via a conic representation.

Lemma 1. 

The product of the voltage and current variables in constraint (6) can be represented as a conic equivalent, which can be relaxed as follows:

$$\begin{array}{c}\hfill ∥\begin{array}{c}2P_{ij,h}\\ l_{ij,h}-v_{i,h}\end{array}∥\le l_{ij,h}+v_{i,h},\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\},\end{array}$$

(17)

where $∥x∥$ represents the $l_2$-norm applied to the vector x.

Proof. 

Equation (6) can be squared on both sides, as follows:

$$\begin{array}{c}\hfill P_{ij,h}^2=v_{i,h}^2{i}_{ij,h}^2,\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\};\\ \hfill P_{ij,h}^2=u_{i,h}{l}_{ij,h},\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\}.\end{array}$$

(18)

Now, the right-hand side of Equation (18) can be represented as a hyperbolic constraint, with the form (18):

$$\begin{array}{c}\hfill P_{ij,h}^2=u_{i,h}{l}_{ij,h}=\frac{1}{4}{\left(u_{i,h}+l_{ij,h}\right)}^2-\frac{1}{4}{\left(u_{i,h}-l_{ij,h}\right)}^2,\left\{\begin{array}{c}\forall ij\in \mathcal{L},\\ \forall h\in \mathcal{H}\end{array}\right\}\\ \hfill 4P_{ij,h}^2+{\left(u_{i,h}-l_{ij,h}\right)}^2={\left(u_{i,h}+l_{ij,h}\right)}^2,\left\{\begin{array}{c}\forall ij\in \mathcal{L},\\ \forall h\in \mathcal{H}\end{array}\right\}.\end{array}$$

(19)

Note that Equation (19) defines a circumference with radius $u_{i,h}+l_{ij,h}$, i.e., it can also be represented as (20):

$$\begin{array}{c}\hfill \sqrt{4P_{ij,h}^2+{\left(u_{i,h}-l_{ij,h}\right)}^2}=u_{i,h}+l_{ij,h},\left\{\begin{array}{c}\forall ij\in \mathcal{L},\\ \forall h\in \mathcal{H}\end{array}\right\},\end{array}$$

(20)

where the left-hand side of (20) can be represented through an $l_2$-norm. Finally, if the equality symbol is relaxed through a lower equal symbol, then constraints (17) and (20) are completely equivalent, and the proof is completed [29].    □

Now, considering the possibility of relaxing the power definition in (6) with (17), the optimization model (1)–(12) becomes a second-order cone programming model, with the structure (21):

$$\begin{array}{cc}\hfill \mathrm{Obj}.\mathrm{Func}.:& \hfill \\ \hfill & min{z}_{\mathrm{approx}}=A_{\mathrm{costs}}+C_{kWh}T\sum _{h\in \mathcal{H}}\sum _{ij\in \mathcal{L}}{R}_{ij}{l}_{ij,h}{\Delta }_h.\hfill \\ \hfill \mathrm{Subject}\mathrm{to}:& \hfill \\ \hfill & P_{ij,h}-R_{ij}{l}_{ij,h}\sum _{k:\left(jk\right)\in \mathcal{L}}{P}_{jk,h}=P_{j,h}^d-P_{j,h}^{pv},\left\{\begin{array}{c}\forall j\in \mathcal{N},\\ j\ne slack,\\ \forall h\in \mathcal{H}\end{array}\right\},\hfill \\ \hfill & u_{j,h}=u_{i,h}-2R_{ij}{P}_{ij,h}+R_{ij}^2{i}_{ij,h},\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\},\hfill \\ \hfill & ∥\begin{array}{c}2P_{ij}\\ l_{ij,h}-v_{i,h}\end{array}∥\le l_{ij,h}+v_{i,h},\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\},\hfill \\ \hfill & z_j{P}_{pv}^{min}\le P_j^{pv}\le z_j{P}_{pv}^{max},\left\{\forall j\in \mathcal{N}\right\},\hfill \\ \hfill & 0\le P_{j,h}^{pv}\le G_h^{pv}{P}_j^{pv},\left\{\forall j\in \mathcal{N},\forall h\in \mathcal{H}\right\},\hfill \\ \hfill & -{\left(i_{ij}^{min}\right)}^2\le l_{ij,h}\le {\left(i_{ij}^{max}\right)}^2,\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\},\hfill \\ \hfill & \sum _{j\in \mathcal{N}}{z}_j\le N_{pv}^{max},\hfill \\ \hfill & z_j\in (0,1),\left\{\forall i\in \mathcal{N}\right\}.\hfill \end{array}$$

(21)

Remark 4. 

Note that, in the optimization model (21), the value of the binary variables $z_j$ was known, as these were obtained after solving the mixed-integer linear programming model (16). In addition, the objective function maintained the component associated with the expected costs of energy losses, given that this factor was essential to ensuring an adequate convergence of the second-order cone programming model (21), as demonstrated by the authors of [24].

3.3. Summary of the Solution Methodology

To summarize the main aspects of the solution methodology, Algorithm 1 presents all of the steps required to solve the problem regarding the optimal placement and sizing of PV generation units in monopolar DC networks via the proposed two-stage solution approach.    

Algorithm 1: General implementation of the proposed two-stage optimization approach.

Data: Obtain the monopolar DC network parameters:

1. Define the slack voltage magnitude as $v_{i,h}=1.0$ pu;

2. Elaborate the optimization model (15) in the CVX programming environment of MATLAB;

3. Solve the mixed-integer linear programming model (15), using the Gurobi solver;

4. Extract the set of nodes where the PV generation units must be installed, i.e., the values of $z_j$;

5. Elaborate the optimization model (21) in the CVX programming environment of MATLAB;

6. Set the values of the binary variables $z_j$ as inputs for the optimization model (21);

7. Solve the second-order cone programming model (21), using the Gurobi solver.

Result: Report the nodes and sizes assigned for the PV generation units.

It is important to note that, in light of the numerical results obtained after implementing the two-stage optimization approach in Algorithm 1, the final size of the PV generation units was reviewed and improved, using the General Algebraic Modeling System (GAMS) software [35].

4. Test Feeder Characteristics

To validate the proposed two-stage optimization methodology, the monopolar DC version of the IEEE 33-bus grid was used as a test feeder. This is a radial distribution network composed of 33 nodes and 32 distribution branches, which is operated with a voltage magnitude of 12,660 V at the terminals of the substation bus. The schematic connections between the nodes of this test feeder are presented in Figure 1.

The parametric information of this test feeder is listed in

Table 1. Parametric information of the IEEE 33-bus grid for monopolar DC studies.

Node i	Node j	${\mathit{R}}_{\mathbf{ij}}$ ($\mathbf{\Omega }$)	${\mathit{P}}_{\mathit{j}}$ (kW)	${\mathit{I}}_{\mathbf{ij}}$ (A)	Node i	Node j	${\mathit{R}}_{\mathbf{ij}}$ ($\mathbf{\Omega }$)	${\mathit{P}}_{\mathit{j}}$ (kW)	${\mathit{I}}_{\mathbf{ij}}$ (A)
1	2	0.0922	100	320	17	18	0.7320	90	20
2	3	0.4930	90	280	2	19	0.1640	90	30
3	4	0.3660	120	195	19	20	1.5042	90	25
4	5	0.3811	60	195	20	21	0.4095	90	20
5	6	0.8190	60	195	21	22	0.7089	90	20
6	7	0.1872	200	95	3	23	0.4512	90	85
7	8	1.7114	200	85	23	24	0.8980	420	70
8	9	1.0300	60	70	24	25	0.8960	420	40
9	10	1.0400	60	55	6	26	0.2030	60	85
10	11	0.1966	45	55	26	27	0.2842	60	85
11	12	0.3744	60	55	27	28	1.0590	60	70
12	13	1.4680	60	40	28	29	0.8042	120	70
13	14	0.5416	120	40	29	30	0.5075	200	55
14	15	0.5910	60	25	30	31	0.9744	150	40
15	16	0.7463	60	20	31	32	0.3105	210	25
16	17	1.2860	60	20	32	33	0.3410	60	20

.

To parameterize the calculation of the objective function (see Equations (2) and (3)) the information reported in

Table 2. Constant parameters regarding the objective function’s parameterization.

Parameter	Value	Unit	Parameter	Value	Unit
$C_{\mathrm{kWh}}$	0.1390	US$/kWh	T	365	days
$t_a$	10	%	$t_e$	2	%
y	20	years	$\Delta h$	1	h
$C_{\mathrm{PV}}$	1036.49	US$/kWp	$C_{\mathrm{O}\&\mathrm{M}}$	0.0019	US$/kWh
$P_i^{pv,max}$	2400	kW	$p_i^{pv,min}$	0	kW
$N_{pv}^{ava}$	3	–	$\Delta V$	±10	%

was employed.

A typical curve for the metropolitan area of Medellín (Colombia) was considered, to evaluate the average daily behavior of the PV generation units. Figure 2 depicts this expected per-unit generation curve. Information regarding solar radiation, irradiance, and temperature were collected from a year’s worth of analysis in the NASA database. This information was processed, with the purpose of obtaining the average behavior for the metropolitan area of Medellín, considering the efficient use of PV modules, in order to obtain the expected generation profile reported in Figure 2. For more details regarding the characterization of these demand and solar curves, please refer to [36].

5. Numerical Results and Discussion

The numerical implementation of the proposed two-stage methodology was carried out in the MATLAB programming environment (version R2021b), running the 64-bit version of Microsoft Windows. Hardware-wise, an Intel(R) Core(TM) i7–7700HQ CPU 2.80 GHz processor and 24 GB of RAM were used. To solve the proposed mixed-integer linear programming and the second-order cone programming models, the CVX convex disciplined tool of MATLAB was used, in addition to its Gurobi solver.

Two simulation scenarios were considered, to validate the effectiveness of the proposed two-stage optimization model:

i.

No thermal limitations in all of the distribution branches in the monopolar DC version of the IEEE 33-bus grid [16];

ii.

The inclusion of the current limitations associated with the conductor sizes assigned to the IEEE 33-bus grid, as defined by the authors of [30].

5.1. Results for the IEEE 33-Bus Grid without Current Limitations

In this subsection, the numerical results obtained by the proposed two-stage approach are compared to three combinatorial optimization methods available in the current literature [16,37]. These correspond to the discrete-continuous versions of the Chu and Beasley genetic algorithm (DCCBGA), the vortex search algorithm (DCVSA), and the generalized normal distribution optimizer (DCGNDO). For all these combinatorial optimizers, a population size of 10 individuals, 1000 iterations, and 100 repetitions were implemented. Note that the proposed method is labeled as a two-stage convex optimizer (TSCO).

Table 3. Evaluation of the comparison methods and the proposed two-stage optimization approach in the monopolar DC version of the IEEE 33-bus grid without thermal constraints.

Method	Site (Node)/Size (kW)	${\mathit{A}}_{\mathbf{cost}}$ (USD/Year)
Bench. case	–	3,644,043.01
CBGA	$\left\{\begin{array}{c}11\left(1162.95\right),14\left(943.48\right),31\left(1482.75\right)\end{array}\right\}$	2,662,724.82
DCVSA	$\left\{\begin{array}{c}9\left(580.31\right),15\left(1291.37\right),31\left(1715.59\right)\end{array}\right\}$	2,662,425.32
GNDO	$\left\{\begin{array}{c}10\left(974.26\right),16\left(920.22\right),31\left(1692.51\right)\end{array}\right\}$	2,662,371.59
TSCO (fixed curve)	$\left\{\begin{array}{c}10\left(974.26\right),16\left(920.22\right),31\left(1692.51\right)\end{array}\right\}$	2,662,371.59
TSCO (variable curve)	$\left\{\begin{array}{c}12\left(1571.78\right),24\left(1676.86\right),30\left(1745.29\right)\end{array}\right\}$	2,561,788.19

presents the numerical results for this simulation scenario, considering two possible operation cases for the PV generation units. The first case considers the operation of the PV plants with maximum power point tracking, as proposed by the authors of [16]. The second approach considers the possibility of generating power without following the expected generation curve, i.e., the PV plants are optimally dispatched in each period. Note that the main idea of these simulation scenarios is to define the best strategy for installing and operating PV plants in electrical networks as a function of the energy resources available, and the dispatch strategy applicable to PV systems [13].

The results in Table 3 show that:

i.

The proposed TSCO approach finds the same numerical solution as the GNDO approach, which confirms that, for the DC version of the IEEE 33-bus grid without thermal limitations in the distribution lines, the best set of nodes to locate PV sources are 10, 16, and 31, with sizes of 974.26, 920.22, and 1692 kW, respectively. This is the best solution for the simulation case where the PV plants follow the maximum point power tracking curve, i.e., they generate the total power available in their terminals though implementing an efficient operation technique [18];

ii.

The DCCBGA and DCVSA are stuck in locally optimal solutions with respect to the DCGNDO and the proposed TSCO approach. The additional gains reached by these methods are USD $353.23$/year and USD $53.73$/year, respectively. However, by comparing the TSCO approach against the benchmark case (operation of the PV network without installed PV generation), it can be observed that the expected reduction in the grid operating costs is USD $\mathrm{981,671.42}$/year, i.e., a significant reduction in the energy investments of the distribution company, with the main advantage that the inclusion of renewables also reduces the carbon footprint.

iii.

The most important finding in Table 3 lies in the comparison between the case involving maximum power point tracking and the optimal dispatch operation scenario. Note that the difference between both cases is about USD $\mathrm{100,583.4}$/year, which corresponds to an additional gain that can be obtained if each PV source is efficiently dispatched in each period. The average energy generation of the PV plants with a fixed curve is about $\mathrm{23,633.04}$ kWh/day, whereas, in the variable generation scenario, $\mathrm{27,976.143}$ kWh/day are produced on average. This result confirms that the efficient operation of PV plants allows for better exploitation of the renewable generation resources available, in comparison to the maximum power point tracking scenario.

To illustrate the positive effect of the variable generation curve on the PV plants, Figure 3 presents a comparative generation curve for the substation bus, regarding the benchmark case and the solutions found by the TSCO approach.

The behavior of the slack generation in Figure 3 shows that: (i) in the benchmark case, the generation output in the slack source follows the aggregated demand generation curve at its terminals (see Figure 2), in conjunction with the total grid power losses; (ii) when the PV generation plants are installed while considering the fixed generation curve scenario, it is noted that, for the period of highest power generation (period 14), the slack generation is zero, unlike the remaining generation periods; and (iii), when the efficient dispatch of PV plants is considered, from period 11 to period 16, the slack generation is zero, which implies that it is indeed possible to maximize the benefit of renewables in monopolar DC networks by this approach. This significant reduction in slack generation is obtained when the PV sources are dispatched without considering the maximum power point, and it shows that, in the studied problem, extracting the total available energy is not always the most economical solution, as the largest sizes regarding renewables allow for better support of power in several periods, which does not occur with maximum power tracking operation.

Note that, in the benchmark case, the slack source generates about $61.5561$ MWh/day. In the maximum power point tracking scenario, this generation is reduced to $37.3198$ MWh/day, whereas, for the variable generation case, the slack source reports a value of about $32.8487$ MWh/day. In addition, it is important to observe that the generation variations in the slack source appear only in the periods with an influence of the PV generation units (periods 7 to 20); in the remaining periods, all three curves are overlapped, as the only available generation source is the slack node. This situation is a clear opportunity for future works to explore the possibility of simultaneously installing energy storage systems to maximize the utilization of renewable generation.

5.2. Results for the IEEE 33-Bus Grid with Current Limitations

This simulation scenario evaluates the effect of considering the current limitations of all the distribution branches on the node location and size selection of PV plants. To demonstrate the effectiveness of the proposed TSCO approach, the exact MINLP model was solved in the GAMS software, and the results were compared for the fixed and variable generation outputs, i.e., the maximum power point tracking scenario and the efficient dispatch of PV sources per period of analysis.

Table 4. Numerical results for the IEEE 33-bus grid while considering current limitations.

Method	Site (Node)/Size (kW)	${\mathit{A}}_{\mathbf{cos}\mathbf{t}}$ (USD/Year)
Bench. case	–	3,644,043.01
GAMS (fixed curve)	$\left\{\begin{array}{c}8\left(1986.84\right),30\left(391.67\right),31\left(939.54\right)\end{array}\right\}$	2,726,761.65
GAMS (variable curve)	$\left\{\begin{array}{c}18\left(482.58\right),25\left(1293.96\right),28\left(2363.82\right)\end{array}\right\}$	2,739,382.36
TSCO (fixed curve)	$\left\{\begin{array}{c}10\left(1327.10\right),24\left(868.21\right),30\left(1326.67\right)\end{array}\right\}$	2,664,816.29
TSCO (variable curve)	$\left\{\begin{array}{c}12\left(1571.78\right),24\left(1676.86\right),30\left(1745.29\right)\end{array}\right\}$	2,561,788.19

compares the GAMS and TSCO results. Note that, in the GAMS software, the BONMIN solver was used to solve the exact MINLP Model (1)–(12) [35].

The results in Table 4 show that:

i.

The implementation of the MINLP model via the BONMIN solver and the GAMS software evidenced that, due to the non-convexities of this optimization model, the solutions reached were only local optima. In the case of the fixed generation curve, the expected reduction concerning the benchmark case was about $25.1721\%$, and, in the case of the variable curve, the reduction only reached a value of $24.8257\%$. Note that these results are counter-intuitive, as a better numerical performance was expected in the case of a variable generation curve; however, this can be explained by the nonlinearities of the exact optimization problem, which causes solvers like BONMIN to be stuck in a locally optimal solution.

ii.

As expected, the solutions reached by the TSCO approach show that variable generation improved the benchmark case by about $29.6993\%$, while the fixed generation curve reported an improvement of about $26.8719\%$. However, the main difference between both solutions corresponded to the modification of one of the nodes, i.e., with regard to the fixed generation case, node 10 changed to node 12 in the variable curve case. Note that the modification in the set of nodal locations was associated with the sensitive behavior of the investment costs when PV generation worked with maximum power point tracking, as compared to the variable operation scenario.

Finally, by comparing the numerical results in Table 3 and Table 4 with respect to the proposed TSCO approach, it is observed that the characteristics of the solution space highly condition the numerical solutions, which implies that the addition of constraints or their simplification will define the final solution. In this sense, more research is required, regarding mixed-integer convex optimization tools to deal with problems pertaining to the installation of renewable generation in electrical networks.

6. Conclusions and Future Work

The problem regarding the optimal placement and sizing of PV generation units in monopolar DC distribution networks was addressed in this research by proposing a two-stage optimization methodology. The first stage relaxed the MINLP model into a mixed-integer linear programming one, by approximating the voltage magnitudes and neglecting the current flowing through the distribution branches. With this mixed-integer linear programming model, the set of nodes to locate the PV sources are selected. The second stage reduced the MINLP model to a nonlinear programming model that could be rewritten as a second-order cone programming one, which defined the optimal size of the PV generation units. Note that these solutions were also refined by solving the exact MINLP model in the GAMS software, setting the binary variables regarding the PV locations.

Two versions of the monopolar DC IEEE 33-bus were considered for numerical validations. The first model neglected the effect of the current bounds in the distribution branches, while the second case considered these thermal limitations. The numerical results obtained from these simulation scenarios showed that the optimal solution for both cases differed by about USD $2444.7$/year, which confirms that including the current effect of the distribution branches can affect the final solution of the problem, given the restrictions regarding branch power flow limitations. This reduces the possible benefit of using large-scale PV systems in some promising nodes, due to the thermal limitations associated with the conductors in the vicinity of these nodes.

As for the free operation of PV generators in each period (i.e., operation without following the maximum power point), significant reductions were observed with respect to the maximum power point tracking case. These reductions were higher than USD $100,583.40$ per year of operation. This reduction in the final objective function demonstrated that the best way to select the size and location of PV generation units in distribution networks is a free power injection in the PV plants during each period, which is defined between zero and the nominal power generation curve when compared to the maximum power point tracking scenario. Note that the effective reduction in the objective function value for the efficient operation scenario was obtained because the conventional generator’s time of use was reduced by at least five hours (for the test feeder), which had a direct impact on the total grid operating costs.

In future works, the following studies can be conducted:

i.

Including, in the proposed optimization model, the set of necessary variables to integrate battery energy storage systems. This, in its conventional formulation, corresponds to a group of linear (convex) constraints associated with the time coupling between the stored energy and the power injection/absorption to/from the monopolar DC network;

ii.

A comparative analysis of different convex optimization with mixed-integer variables, to solve the MINLP problem and compare their performances against that of the proposed TSCO approach. Note that some mixed-integer convex formulations using second-order cone programming theory can directly solve the problem regarding the sizing and location of PV plants simultaneously (in one stage), which constitutes an excellent opportunity to validate the effectiveness of the proposed TSCO approach, especially in large-scale monopolar DC networks;

iii.

Considering the demand and PV generation curves, including the uncertainties associated with the stochastic nature of these variables and their effect on the final grid operation plan, as well as some physical constraints associated with the space available for the installation of PV sources, as there are some nodes where this is a limitation;

iv.

The reformulation of the TSCO approach, using convex approximations that allow for simultaneously dealing with radial and meshed distribution networks. These reformulations can be based on the nodal voltage method and semi-definite programming models.