Optimal Reactive Power Compensation via D-STATCOMs in Electrical Distribution Systems by Applying the Generalized Normal Distribution Optimizer

Abstract

This research deals with the problem regarding the optimal siting and sizing of distribution static compensators (D-STATCOMs) via the application of a master–slave optimization technique. The master stage determines the nodes where the D-STATCOMs must be located and their nominal rates by applying the generalized normal distribution optimizer (GNDO) with a discrete–continuous codification. In the slave stage, the successive approximations power flow method is implemented in order to establish the technical feasibility of the solution provided by the master stage, i.e., voltage regulation and device capabilities, among other features. The main goal of the proposed master–slave optimizer is to minimize the expected annual operating costs of the distribution grid, which includes the energy loss and investment costs of the D-STATCOMs. With the purpose of improving the effectiveness of reactive power compensation during the daily operation of the distribution grid, an optimal reactive power flow (ORPF) approach is used that considers the nodes where D-STATCOMs are located as inputs in order to obtain their daily expected dynamical behavior with regard to reactive power injection to obtain additional net profits. The GNDO approach and the power flow method are implemented in the MATLAB programming environment, and the ORPF approach is implemented in the GAMS software using a test feeder comprising 33 nodes with both radial and meshed configurations. A complete comparative analysis with the Salp Swarm Algorithm is presented in order to demonstrate the effectiveness of the proposed two-stage optimization approach in the fixed operation scenario regarding the final objective function values. In addition, different tests considering the possibility of hourly power injection using D-STATCOMs through the ORPF solution demonstrate that additional gains can be obtained in the expected annual operative costs of the grid.

1. Introduction

Medium- and low-voltage distribution networks have undergone important changes in their physical and computational structures that have transformed these passive grids into to active distribution networks [1,2]. These changes have been provoked by new regulations regarding the integration of renewable energy resources that have focused on the decabornization of the world’s energy matrix [3,4,5]. The main goal of electricity distribution networks is to transport electricity to end-users with the necessary standards of efficiency, quality, and reliability. This requires minimizing energy losses and improving transport processes in order to achieve more efficient and competitive networks [6].

To improve the efficiency of electrical distribution networks with regard to reducing power losses, reactive power compensation is one of the most widely recognized methods, given its contribution in this regard and other benefits such as power factor correction, increased transport capacity, improved operation of network lines and devices, voltage stability, and improved voltage profiles. All of these depend on different operating constraints [7], as it is of great importance to reduce losses and minimize investment costs for current distribution systems.

One of the best options for energy loss reduction is distributed generation (DG). The authors of [8] presented a decision-making analysis aimed at determining the optimal location and sizes for DG in radial distribution networks, considering the improvement of the voltage profiles and the minimization of the total grid’s active and reactive power losses. However, DG may have very high initial installation costs, especially when compared to strategies such as grid reconfiguration and shunt reactive power compensation with capacitors [9]. Capacitive compensation can yield the benefits of loss reduction, power factor correction, and voltage profile improvement to the fullest if the location and size are efficiently determined. However, capacitor banks typically inject reactive power in fixed steps [10], which reduces the positive effect of reactive power injection on reducing energy losses, mainly when daily demand profiles with high variations are considered in the grid operation environment.

In order to improve the advantages of using reactive power to minimize total grid power losses, it is possible to use distribution static compensators (D-STATCOMs). These are regulation devices based on a power electronics voltage source converter that add flexibility to the power distribution network. In comparison with a variable-step capacitor bank, a D-STATCOM is a more flexible device, as it injects the exact amount of reactive power in accordance with the network requirements. This is by means of efficient control techniques [11]. It is important to mention that D-STATCOMs can inject and absorb reactive power with very fast dynamic responses through the injection of phase-shifted current to the system at the common coupling point, which aids in power factor correction and harmonic filtering, among other benefits [12]. D-STATCOMs can be located in any of the nodes of the distribution system, so it is important to determine their optimal siting and sizing using efficient optimization techniques for the reduction of energy losses and the minimization of operating costs.

In the specialized literature, there are multiple reports regarding the optimal location and sizing of D-STATCOMs in distribution networks. Some of these approaches are discussed in this document. The authors of [12] presented a complete way to determine the optimal location and size of D-STATCOMs via analytical and heuristic optimization methods. In addition, they presented the typical objective functions of the specialized literature to improve network performance, i.e., voltage stability and power losses indicators. The work by [13] proposed a multi-objective particle swarm optimizer to place and size D-STATCOMs while considering simultaneous power grid reconfiguration. As objective functions, the minimization of active power losses, the voltage stability index, and the load capacity factor of the distribution lines were considered. The main feature of this approach is that the optimization process was carried out only under maximum load conditions, which is not considered a suitable scenario given the potential oversizing of the compensating devices, since the consumption of active and reactive power are variable inputs.

The authors of [14] proposed a multi-objective fuzzy approach based on ant colony optimization in order to solve the simultaneous and proximate reconfigurations (sizing and placement) of PV sources and D-STATCOMs in distribution systems. The objective was to minimize network losses and to improve voltage profiles and feeder load balancing operation characteristics. This methodology was validated in the IEEE 33-bus grid with excellent numerical results. The authors of [15] presented a heuristic method based on voltage and power loss indexes in order to locate and size D-STATCOMs in radial electrical distribution networks. Numerical validation of this heuristic approach was performed on the IEEE 33-bus test feeder. However, the authors only considered maximum load conditions.

Table 1. Main literature reports regarding the optimal placement and sizing of D-STATCOMs in distribution grids.

Solution Method	Objective Function	Ref.	Year
Evolutionary programming	Power loss minimization and voltage profile improvement	[16]	2011
Particle swarm optimization	Power loss minimization and voltage profile improvement	[17]	2014
Sensitivity index	Power loss minimization and voltage profile improvement	[18]	2015
Gravitational search algorithm	Power loss minimization, voltage profile improvement, and investment cost minimization	[19]	2016
Adaptive particle swarm optimization	Energy loss and investment cost minimization	[20]	2018
Particle swarm optimization	Energy loss minimization	[21]	2018
Salp swarm algorithm	Energy loss and investment cost minimization	[22]	2022
Two-stage convex optimization model	Energy loss and investment cost minimization	[23]	2022

summarizes the different algorithms used in the literature to solve the problem concerning the placement and sizing of D-STATCOMs in distribution networks.

The main features of the optimization methods in Table 1 are the following: (i) the most common function is the minimization of energy losses, and (ii) the studied problems show two tendencies: the first is related to metaheuristics, and the second offers convex formulations or approaches that combine convex and combinatorial methods.

Based on the aforementioned literature review, the problem of reactive power compensation in power systems is still an area of interest for both academia and the industry. Therefore, this research article proposes the following contributions:

• We provide a new solution method based on application of the generalized normal distribution optimizer (GNDO) for locating and sizing D-STATCOMs in distribution networks with radial and meshed topologies while using a discrete–continuous codification.
• We combine the GNDO approach with an efficient power flow multi-period approach that allows solving of the technical constraints of the optimization problem, i.e., power balance, voltage regulation, and device capabilities, among others.
• We improve the final solution obtained with the proposed master–slave optimizer by using the set of nodes where the D-STATCOMs must be located as inputs for the optimal reactive power flow (ORPF) problem, aiming to further minimize the final expected annual operating costs of the distribution grid.

Note that the selection of the GNDO as a solution technique to address the optimal location and sizing problem for D-STATCOMs in distribution networks is based on three facts:

• The GNDO is a metaheuristic method inspired by the classical theory of normal probability distributions. It considers an initial population that evolves throughout the iterative process, considering the means and the standard deviation as advance parameters.
• The computational implementation of the GNDO is simple and requires only a few mathematical programming skills to adapt it to any optimization problem that includes binary and continuous variables.
• Multiple reports in the specialized literature have confirmed that the GNDO approach is efficient at solving complex optimization problems such as the placement location of renewable energy resources in AC and DC networks [24,25] and parameter extraction for photovoltaic models [26], among other optimization problems.

Note that the proposed solution methodology addressed in this research to locate, size, and operate D-STATCOMs in electrical distribution networks is different from recently published literature reports, since the GNDO approach allows reaching better objective function values with respect to the literature in the case of fixed reactive power injection, and also, with the improvement stage based on the ORPF solution, we find additional profits for the utility company that have not been previously presented by solution methodologies based on combinatorial optimizers available in the current literature. It is important to mention that this work focuses on the optimal design and operation of D-STATCOMS in distribution networks under normal operating conditions and upon the basis of the demand curves provided by utilities. However, as future research, an interesting development could consider the probabilistic nature of the loads and new paradigms regarding the transformation of energy consumption habits with the massive integration of electric vehicles into medium- and low-voltage distribution networks [27].

The remainder of this work has the following structure presented: Section 2 presents the mathematical formulation of the D-STATCOM optimal integration problem in power distribution systems via mixed-integer nonlinear programming models. Section 3 presents the proposed solution method, which is based on the GNDO and the successive approximations power flow approach, with the aim of determining the nodes where the D-STATCOMs must be located as well as their sizes. These nodes are also used as inputs for the ORPF formulation to reach additional reductions in the objective function value. Section 4 outlines the main characteristics of the test system, which is a 33-node IEEE system that includes demand scenarios for residential, industrial, and commercial users. Section 5 presents the results of the simulation as well as its analysis and discussion. Finally, Section 6 describes the main conclusions derived from this study.

2. Mathematical Formulation

This section formulates the problem regarding the optimal location and size of D-STATCOMs in electrical distribution networks as a general mixed-integer nonlinear programming (MINLP) model. The continuous variables of the model denote the voltages, angles, and active and reactive power variables, while the discrete variables correspond to the possibility of integrating D-STATCOMs in all the nodes of the distribution network. In addition, these equations take a nonlinear form due to the products between them and the trigonometric functions in the power balance equations. The complete MINLP model of the studied problem is presented below.

2.1. Objective Function Structure

The problem under study is addressed in this research from an economic perspective. The main idea of installing these devices is to minimize the expected operating cost of the network, which is mainly associated with the energy loss costs for a one-year period of operation ($f_1$), to which the annualized investments in D-STATCOMs ($f_2$) are added.

$$\begin{array}{c}{f}_1=C_{kWh}T\sum _{h\in \mathcal{H}}\sum _{k\in \mathcal{N}}\sum _{m\in \mathcal{N}}{Y}_{km}{V}_{kh}{V}_{mh}cos\left(\begin{array}{c}{\delta }_{kh}-{\delta }_{mh}-{\theta }_{km}\end{array}\right){\Delta }_h,\end{array}$$

(1)

$$\begin{array}{c}{f}_2=T{\displaystyle \frac{k_1}{k_2}}\sum _{k\in \mathcal{N}}\left(\alpha {\left(Q_k^{st}\right)}^2+\beta Q_k^{st}+\gamma \right)Q_k^{st},\end{array}$$

(2)

where for $f_1$ and $f_2$, $C_{kWh}$ corresponds to the average cost of energy in $kW/h$; T represents a constant associated with the number of days in a regular year (length of the period of study); $Y_{km}$ defines the magnitude of the admittance matrix that relates k node to m node, which also has an associated angle ${\theta }_{km}$; $V_{kh}$ and $V_{mh}$ are the voltages that relate nodes k and m in period of time h with voltage angles defined as ${\delta }_{kh}$ and ${\delta }_{mh}$, respectively; ${\Delta }_h$ is a period of time assigned as $1h$; $\mathcal{H}$ and $\mathcal{N}$ refer to the set of time periods and nodes in the network, respectively; $k_1$ and $k_2$ are two positive parameters that denote the annualization costs of the D-STATCOMs (i.e., $k_1$ represents the annual investment required by the D-STATCOM, and $k_2$ refers to the useful life of the device); and the parameters $\alpha$, $\beta$, and $\gamma$ are positive scalars that represent the installation costs of the D-STATCOMs with nominal power and generation capacity $Q_k^{st}$.

It is worth mentioning that Equation (1) is the most typical representation of power and energy losses in distribution networks, and it can be expressed as the difference between the energy input (energy injected by power sources) and energy output (energy consumed by constant power loads) [23]. However, this is typically represented by the product between voltages and trigonometric functions of their angles, as this allows evidencing the effect of the power balance constraints (which are presented later in this document) regarding variations in active or reactive power injection. In the case of the quantification of the investment costs of D-STATCOMs, as defined by Equation (2), this function was adapted from [19], where the authors presented different cubic functions to represent the costs of variable reactive power compensators based on power electronics for power system applications.

The objective function corresponding to the object of minimization is defined in this research as the sum of the components $f_1$ and $f_2$, as presented in Equation (3) for the variable $A_f$.

$$\begin{array}{c}min{A}_f=f_1+f_2,\end{array}$$

(3)

where $A_f$ is denoted as the objective function value to be minimized.

Remark 1. 

The structure of the objective function (3) is nonlinear and non-convex due to the following facts: (i) the component $f_1$ is defined as a product between the continuous variables (voltages) and the trigonometric function of the voltage angles, and (ii) the component $f_2$ is a cubic function defined with the expected sizes of the D-STATCOMs, which makes the final function value $A_f$ a complex nonlinear objective function that requires advanced optimization techniques in order to obtain a solution.

In this research, to solve the complex objective function regarding the installation and sizing of D-STATCOMs in distribution grids with radial or meshed structures, the GNDO approach is combined with the successive approximations power flow method. In addition, an improving stage based on the ORPF problem is implemented in order to obtain additional reductions in the final $A_f$ by fixing the location of the D-STATCOMs determined via the GNDO approach.

2.2. Set of Constraints

To establish the location and size of D-STATCOMs in radial and meshed distribution networks, it is necessary to observe several constraints regarding the physical operation of the distribution network. Some of these have to do with active and reactive power balance, voltage grid regulation characteristics, and the sizes of the reactive power compensation devices. These constraints are listed below and defined for $\forall k\in \mathcal{N}$ and $\forall h\in \mathcal{H}$.

$$\begin{array}{c}{P}_{kh}^g-P_{kh}^d=\sum _{m\in \mathcal{N}}{Y}_{km}{V}_{kh}{V}_{mh}cos\left({\delta }_{kh}-{\delta }_{mh}-{\theta }_{km}\right),\end{array}$$

(4)

$$\begin{array}{c}{Q}_{kh}^g-Q_{kh}^d+Q_k^{st}=\sum _{m\in \mathcal{N}}{Y}_{km}{V}_{kh}{V}_{mh}sin\left({\delta }_{kh}-{\delta }_{mh}-{\theta }_{km}\right),\end{array}$$

(5)

$$\begin{array}{c}{V}_{min}\le V_{k,h}\le V_{max},\end{array}$$

(6)

$$\begin{array}{c}{Z}_k{Q}_{min}^{st}\le Q_k^{st}\le Z_k{Q}_{max}^{st},\end{array}$$

(7)

$$\begin{array}{c}\sum _{k\in \mathcal{N}}{Z}_k\le N_A^{st},\end{array}$$

(8)

where $Q_k^{st}$ is the reactive injection generated by the D-STATCOM at node k; the variables $P_{kh}^g$ and $Q_{kh}^g$ are the active and reactive power injection in power sources connected to k in the period h; the parameters $Q_{kh}^d$ and $P_{kh}^d$ are the dynamical loads associated with node k at time h, which are defined as a function of the residential, industrial, and commercial users connected to the distribution grid; $Z_k$ is a binary variable (i.e., it takes values of 0 or 1) that indicates the connection or not of a D-STATCOM at node k; and $N_A^{st}$ corresponds to the maximum number of D-STATCOMs that can be installed in the distribution network.

Constraints (4)–(7) can be interpreted as follows. Equations (4) and (5) defines the power balance constraints at each node of the network for each period. These are known as the power flow constraints [28]. Inequality constraint (6) refers to the voltage regulation imposed by the regulatory entities for the adequate operation of the distribution network [29]. Inequality constraint (7) concerns the upper and lower bounds assigned for the D-STATCOMs that are to be installed in the distribution network. Finally, inequality constraint (8) is associated with the maximum number of D-STATCOMs that can be placed in the distribution grid.

3. Proposed Solution Method

To solve the complex MINLP formulation in Equations (1)–(8), this research presents a two-stage optimization model. In the first stage, a master–slave optimization algorithm based on the GNDO approach and the successive approximations power flow method is implemented in order to determine the nodes where the D-STATCOMs must be placed as well as their initial sizing. In the second stage, an approach of improvement is implemented, which involves solving the ORPF problem by fixing the nodes where the D-STATCOMs must be located and recalculating their sizes. Each of the previous stages and algorithms are presented in this section.

3.1. Master–Slave Optimization Algorithm

This subsection discusses the master–slave approach, which comprises the GNDO and the successive approximations power flow method.

Slave Stage: Successive Approximations Power Flow Method

This power flow method is a well-known approach that reformulates Equations (4) and (5) using a complex equivalent in order to find a solution for radial and meshed distribution networks. This method was initially proposed by [30], who were inspired by the classical Gauss–Seidel power flow approach. The general iterative power flow formula for the successive approximations power flow method is defined in Equation (9). Note that the dimensions of vectors and matrices for the successive approximations power flow method and a numerical example can be consulted in [31].

$$\begin{array}{c}{\mathbb{V}}_{dh}^{m+1}={\mathbb{Y}}_{dd}^{-1}\left[{\mathrm{diag}}^{-1}\left({\mathbb{V}}_{dh}^{m,★}\right)\left({\mathbb{S}}_{sth}^{★}-{\mathbb{S}}_{dh}^{★}\right)-{\mathbb{Y}}_{ds}{\mathbb{V}}_{sh}\right],\left\{h\in \mathcal{H}\right\}\end{array}$$

(9)

where ${\mathbb{V}}_{dh}$ is the vector containing all the voltage demand variables for the period of time h; ${\mathbb{V}}_{sh}$ is the vector that contains the slack voltage output in the period of time h; ${\mathbb{Y}}_{dd}$ represents a square matrix containing the admittance relation of the demand nodes; ${\mathbb{Y}}_{ds}$ is a rectangular matrix that associates the demand nodes with the slack source regarding their admittance interconnections; ${\mathbb{S}}_{sth}^{★}$ is a vector that contains all the reactive power injections provided by the D-STATCOMs in the nodes where they are connected (note that the reactive power injection is zero at nodes where there are no D-STATCOMs); and ${\mathbb{S}}_{dh}^{★}$ is a complex vector comprising the constant power consumption of the distribution network per period of time h. Note that m represents the iteration counter, and the recursive power flow formula in Equation (9) converges if the following criterion is met:

$$max\left\{\left|\left|{\mathbb{V}}_{dh}^{m+1}\right|-\left|{\mathbb{V}}_{dh}^m\right|\right|\right\}\le ϵ.$$

(10)

It is worth mentioning that in Equation (10), the stopping criterion has to be reached before the maximum number of iterations. Here, $ϵ$ is the assigned tolerance.

Remark 2. 

The main advantage of using the successive approximations power flow method to deal with the power flow problem in radial and meshed distribution grids is that it ensures convergence to the solution if the grid is connected (there are no isolated areas) and is operating far from the voltage collapse point. Its convergence can be demonstrated by applying the Banach fixed-point theorem [30].

Once the recursive power flow formula (9) has converged to the power flow solution, the power generation in the slack source can be calculated as presented in Equation (11).

$${\mathbb{S}}_{sh}^{*}=\mathrm{diag}\left({\mathbb{V}}_{sh}^{*}\right)\left[{\mathbb{Y}}_{ss}{\mathbb{V}}_{sh}+{\mathbb{Y}}_{sd}{\mathbb{V}}_{dh}\right]\forall \left\{h\in \mathcal{H}\right\}.$$

(11)

With this solution, it is possible to evaluate the expected annual energy loss cost of the distribution grid ($f_1$), as stated in the objective function. In addition, the constraint regarding voltage regulation is checked, (i.e., see inequality constraint (6)). If it is fulfilled, then a zero value is added to the objective function; if not, then a penalization regarding the amplitude of the voltage deviation with respect to its upper and lower bounds is added to the objective function value. Note that adding penalty factors to the objective function is a typical approach employed to find feasible solutions while exploring and exploiting the solution space. More details about penalty factors in metaheuristics can be consulted in [32].

3.2. Master Stage: GNDO

The GNDO approach is a recently developed metaheuristic optimization approach based on the behavior of the Gaussian distribution that can model multiple physical and natural phenomena with a high level of precision [26]. This metaheuristic optimization method is part of the mathematics-inspired approaches, and it uses two main parameters to explore and exploit the solution space: the mean value ($\mu$) and the standard deviations ($\delta$). In addition, the GNDO approach considers two main steps during this process: local and global searching, which are discussed below.

3.2.1. Local Exploration

Local exploration is a common search technique in metaheuristic-based optimization methods, the main idea of which is to explore in the vicinity of a current solution (i.e., minimum variations) in order to find a better one. In other words, it is an exploitation stage. The rule defined in the GNDO for this exploration is defined below:

$$\begin{array}{c}\hfill v_i^t={\mu }_i+{\delta }_i\eta ,i=1,2,\cdots ,N_i,\end{array}$$

(12)

where $v_i^t$ corresponds to the trailing vector associated with the current solution i at iteration t; ${\mu }_i$ denotes the generalized mean location of the individual solution i during iteration t; ${\delta }_i$ corresponds to the generalized standard deviation of this solution; and $\eta$ is a penalty factor. It is worth mentioning that $N_i$ is a constant parameter related to the number of individuals that make up the algorithm population.

For the optimal location and sizing problem regarding D-STATCOMs in radial and meshed distribution networks, the proposed codification is a discrete–continuous representation, where the integer part defines the nodes where the compensators must be installed, and the continuous part of the codification corresponds to their optimal sizes. Equation (13) presents the proposed codification.

$$\begin{array}{c}\hfill x_i^t=\left[15,k,\cdots ,28|1.3566,Q_k^{st},\cdots ,2.0101\right]\end{array}$$

(13)

Remark 3. 

It is important to observe that with the codification presented in Equation (13), the objective function component related to the D-STATCOM investment costs ($f_2$) is completely known.

To determine the value of the parameters ${\mu }_i$, ${\delta }_i$, and $\eta$ defined in Equation (12), the following set of equations is used:

$$\begin{array}{c}\hfill {\mu }_i=\frac{1}{3}\left(x_i^t+x_{\mathrm{best}}^t+M\right),\end{array}$$

(14)

$$\begin{array}{c}\hfill {\delta }_i=\frac{1}{\sqrt{3}}{\left({\left(x_i^t-\mu \right)}^2+{\left(x_{\mathrm{best}}^t-\mu \right)}^2+{\left(M-\mu \right)}^2\right)}^{\frac{1}{2}},\end{array}$$

(15)

$$\begin{array}{c}\hfill \eta =\left\{\begin{array}{ccc}{\left(-\log \left({\lambda }_1\right)\right)}^{\frac{1}{2}}cos\left(2\pi {\lambda }_2\right)& \mathrm{if}& a\le b\\ {\left(-\log \left({\lambda }_1\right)\right)}^{\frac{1}{2}}cos\left(2\pi {\lambda }_2+\pi \right)& \mathrm{if}& a>b\end{array}\right.\end{array}$$

(16)

In Equations (14)–(16), it can be noted that:

• The parameters a, b, ${\lambda }_1$, and ${\lambda }_2$ are values generated with a uniform distribution in the interval $\left[01\right]$;
• $x_{\mathrm{best}}^t$ represents the best current solution reached as of iteration t;
• M is a vector that represents the average (mean value) associated with all the individuals in the current population.

To calculate the vector M, Equation (17) is used.

$$\begin{array}{c}\hfill M=\frac{1}{N_i}\sum _{i=1}^{N_i}{x}_i^t.\end{array}$$

(17)

3.2.2. Global Exploration

Global exploration is an intrinsic characteristic of combinatorial optimization methods, where the main idea is to explore the solution space in search of promissory solution regions that may contain the global optimum of the studied problem [24]. The global exploration stage, as proposed by the authors of [26], is defined in Equation (18).

$$\begin{array}{c}\hfill v_i^t=x_i^t+\beta \times \left(|{\lambda }_3|\times v_1\right)+\left(1-\beta \right)\times \left(|{\lambda }_4|\times v_2\right),\end{array}$$

(18)

where $\beta \times \left(|{\lambda }_3|\times v_1\right)$ is a component associated with information regarding the local exploration, and $\left(1-\beta \right)\times \left(|{\lambda }_4|\times v_2\right)$ is a component related to the global exploration stage. Note that ${\lambda }_3$ and ${\lambda }_4$ denote random numbers generated with a normal distribution, and $\beta$ is an adjusting parameter contained in the interval $\left[01\right]$, which is randomly obtained using a uniform distribution. In addition, $v_1$ and $v_2$ are also two trail vectors obtained with the following rules:

$$\begin{array}{c}\hfill v_1=\left\{\begin{array}{ccc}{x}_i^t-x_j^t& \mathrm{if}& A_f\left(x_i^t\right)<A_f\left(x_j^t\right)\\ x_j^t-x_i^t& & \mathrm{otherwise}\end{array}\right.\end{array}$$

(19)

$$\begin{array}{c}\hfill v_2=\left\{\begin{array}{ccc}{x}_k^t-x_k^t& \mathrm{if}& A_f\left(x_k^t\right)<A_f\left(x_m^t\right)\\ x_m^t-x_k^t& & \mathrm{otherwise}\end{array}\right.\end{array}$$

(20)

Note that subscripts j, k, and m correspond to integer numbers between 1 and $N_i$ that are associated with three different individuals contained in the population $x^t$. These subscripts are different from each other as well as from the current individual i.

Prior to deciding which $v_i^t$ will be part of the next population of individuals, it is mandatory to check its feasibility regarding the upper and lower limits of the decision variables, i.e.,

$$\begin{array}{c}\hfill v_{i,l}^t=\left\{\begin{array}{ccc}{v}_{i,l}^t& \mathrm{if}& x_l^{min}\le v_{i,l}^t\le x_l^{max}\\ v_{i,l}^t=x_l^{min}+{\lambda }_5\left(x_l^{max}-x_l^{min}\right)& & \mathrm{otherwise}\end{array}\right.\end{array}$$

(21)

where $x_j^{min}$ and $x_j^{max}$ are the lower and upper admissible bounds of variable j, respectively. Note that due to the integer nature of the variables regarding the nodes where the D-STATCOMs must be located, the first $N_A^{st}$ position of each solution vector is rounded to the nearest integer value. Finally, to select the solution that will be part of the next generation of individuals, the following selection rule is applied:

$$\begin{array}{c}\hfill x_i^{t+1}=\left\{\begin{array}{ccc}{v}_i^t& \mathrm{if}& A_f\left(v_i^t\right)<A_f\left(x_i^t\right)\\ x_i^t& & \mathrm{otherwise}\end{array}\right.\end{array}$$

(22)

3.2.3. Implementation of the GNDO Approach

As explained by Equations (9)–(22), to illustrate the application of the GNDO approach to the problem regarding the optimal placement and sizing of D-STATCOMs in electrical distribution grids, Algorithm 1 is used.

Algorithm 1: General application of the GNDO approach to select nodes and define sizes for D-STATCOMs in distribution networks.

Note that Algorithm 1 uses all of the formulas involved in the GNDO approach (9)–(22) in order to explore and exploit the solution space by dividing the search into local and global explorations.

3.3. Optimal Reactive Power Flow Solution

This improvement stage deals with defining the optimal sizes of the D-STATCOMs in the case of daily operation, i.e., considering variable reactive power injection throughout the operation scenario. To this effect, the first $N_A^{st}$ positions of the final best solution $x_{\mathrm{best}}^t$, i.e., the nodes where the D-STATCOMs must be located, are fixed in the optimization model (1)–(8), which transforms the exact MINLP model into an NLP one known as the optimal reactive power flow problem [33].

To summarize the proposed improvement, a flow diagram is presented in Figure 1.

Remark 4. 

Note that the solution of the NLP problem in Figure 1 is independent of the optimization tool. However, this research resorts to the GAMS software, as it has been used in multiple literature reports with excellent numerical performance for power system optimization with low computational effort [34].

4. Test System Characterization

To demonstrate the efficiency and robustness of the proposed optimization, the IEEE 33-bus grid with radial and meshed configurations was employed [22]. The electrical configuration of this distribution grid is depicted in Figure 2. Note that this electrical system operates at substation terminals with a nominal voltage of $12.66$ kV. In addition, the parametric information regarding the peak load consumption and distribution lines is presented in

Table 2. Load and distribution line parameters of the IEEE 33-bus grid.

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

.

To consider the effect of load classification in the IEEE 33-bus grid for optimally locating and sizing the D-STATCOMs, the electrical curves of residential, industrial, and commercial users are shown in

Table 3. Expected percentages of consumption for residential, industrial, and commercial users.

Hour	Ind. (pu)	Res. (pu)	Com. (pu)
1	0.56	0.69	0.20
2	0.54	0.65	0.19
3	0.52	0.62	0.18
4	0.50	0.56	0.18
5	0.55	0.58	0.20
6	0.58	0.61	0.22
7	0.68	0.64	0.25
8	0.80	0.76	0.40
9	0.90	0.90	0.65
10	0.98	0.95	0.86
11	1	0.98	0.90
12	0.94	1	0.92
13	0.95	0.99	0.89
14	0.96	0.99	0.92
15	0.9	1	0.94
16	0.83	0.96	0.96
17	0.78	0.96	1
18	0.72	0.94	0.88
19	0.71	0.93	0.76
20	0.70	0.92	0.73
21	0.69	0.91	0.65
22	0.67	0.88	0.5
23	0.65	0.84	0.28
24	0.60	0.72	0.22

.

To calculate the objective function value regarding component $f_2$ (see Equation (2)), all the parameters listed in

Table 4. Parametric information to evaluate the objective function regarding the investment costs of the D-STATCOMs.

Par.	Value	Unit	Par.	Value	Unit
$C_{\mathrm{kWh}}$	0.1390	USD*kWh	T	365	Days
${\Delta }_h$	0.50	h	$\alpha$	0.30	USD/MVAr${}^3$
$\beta$	−305.10	USD*MVAr${}^2$	$\gamma$	127,380	USD/MVAr
$k_1$	6/2190	1/Days	$k_2$	10	Years

are used [19,35]. Note that the evaluation of this objective function considers that the variable $Q_k^{st}$ is defined in Mvar.

In this research, the concept of a benchmark case is assigned to the scenario where the distribution system is operated without considering shunt compensation, i.e., it is the simulation scenario where the total active and reactive power consumptions are always supplied by the generation source, which, in the case of distribution networks, is the substation bus when there is no penetration of distributed energy resources. In the case of the IEEE 33-bus grid with a radial connection, the benchmark case is calculated with the parameters in Table 4, which corresponds only to an $f_1$ value of about USD $130,613.90$ per year of operation, with $f_2$ being equal to zero since no D-STATCOMs are considered.

5. Results and Discussion

To validate the proposed solution methodology, the solution of the optimization model for locating and sizing D-STATCOMs in electrical distribution via the AESDC was implemented in MATLAB version R2021b on a computer with an Intel Core i7-10750H @2.6 GHz and 16.0 GB of DDR4 RAM at 2300 MHz on a 64-bit version of Microsoft Windows 10 Home. The validation of the proposed master–slave optimization model in conjunction with the improvement stage was carried out while considering the following facts:

• Two configurations are considered for the IEEE 33-bus grid: a radial configuration with the original structure presented in Figure 2 and a meshed configuration with three tie-lines added to the radial topology, as proposed in [22].
• The GNDO approach, combined with the successive approximations power flow method, is used to obtain the initial nodal locations and sizes for the D-STATCOMs with fixed reactive power injection. Thus, the location of these devices is fixed in the exact MINLP model while aiming for additional gains for operating these reactive power compensators with variable reactive power throughout the day.
• For comparison, the MINLP model’s exact solutions obtained with the BONMIN and COUENNE solvers and the recently developed Salp Swarm Algorithm (SSA) [22] are considered in order to confirm the effectiveness and robustness of the proposed optimization approach.

Note that in order to ensure the effectiveness of the proposed master–slave optimization regarding the solution quality, 100 consecutive evaluations of the complete solution methodology were considered, with 1000 iterations for each of them. With these values, a statistical study of the GNDO performance was conducted considering mean, maximum, minimum, and standard deviation parameters. It is worth mentioning that in order to to make fair comparisons with the SSA approach presented by the authors of [22] as well as the implementation of the exact model in GAMS software with the BONMIN and COUNNE solvers, all results were obtained with our own implementations that were conducted during this research.

It is worth mentioning that in order to ensure that all the numerical validations presented in this paper were correct, the following validation procedure was followed.

• The solution of the multi-period optimal power flow problem for the benchmark case was implemented in MATLAB with our own scripts. However, to ensure the effectiveness of this approach, a comparative analysis with the DIgSILENT software was performed.
• The convergence analysis of the power flow problem for the proposed successive approximations power flow method was proven by the authors of [30], which implies that under the studied conditions this method always converges to the power flow solution.
• The solution of the MINLP model with he GNDO approach combined with the successive approximations method was also validated via GAMS and DIgSILENT.

5.1. Radial Configuration Results

Table 5. Optimal solutions obtained by the studied optimization methods for a radial grid configuration.

Method	Nodes	Sizes (Mvar)	Annual Costs (USD/year)	Red. (%)
Ben. Case	—	—	130,613.90	0.00
COUENNE	$\left[\begin{array}{ccc}5& 6& 11\end{array}\right]$	$\left[\begin{array}{ccc}0.0000& 0.53600& 0.27466\end{array}\right]$	115,960,99	11.22
BONMIN	$\left[\begin{array}{ccc}8& 25& 30\end{array}\right]$	$\left[\begin{array}{ccc}0.29796& 0.09204& 0.51265\end{array}\right]$	109,560.85	16.12
SSA	$\left[\begin{array}{ccc}13& 25& 30\end{array}\right]$	$\left[\begin{array}{ccc}0.25851& 0.10547& 0.52801\end{array}\right]$	108,249.36	17.12
GNDO	$\left[\begin{array}{ccc}14& 25& 30\end{array}\right]$	$\left[\begin{array}{ccc}0.23083& 0.09996& 0.53905\end{array}\right]$	108,215.94	17.15

shows the comparison between the SSA approach, the BONMIN and COUENNE solvers, and to the proposed GNDO approach that considers fixed reactive power injection throughout the day for the IEEE 33-bus system with a radial configuration.

The numerical results presented in Table 5 allow stating that:

• The GAMS solvers (COUENNE and BONMIN) are stuck in locally optimal solutions, with reductions of 11.22 and 16.12%, respectively, with respect to the benchmark case. This can be attributed to the fact that the MINLP structure of the model (1)–(8) causes exact solution methods based on Branch and Bound and combined with interior point methods to get stuck in local solutions without the ability to escape from them because of the non-convexity of the solution space.
• The proposed GNDO approach improved the best solution reported in the literature, which was obtained with the SSA approach, i.e., by about USD $33.42$ per year of operation. This improvement was reached since the GNDO detected that node 14 is a better location for one of the D-STATCOMs, instead of node 13 as per the SSA approach.
• The total installed reactive power obtained with the GNDO approach was 869.84 kvar, while the solution found with the SSA approach installed 891.99 kvar. This implies that the proposed GNDO approach improved the objective function value by selecting a better set of nodes to place the D-STATCOMs, with the main advantage being that less investment is required for these devices.

To demonstrate that the hourly variation of the reactive power injection with D-STATCOMs in distribution networks allows additional gains in the expected annual grid operating costs,

Table 6. Additional improvements reached when variable reactive power is used for daily compensation with D-STATCOMs for the radial grid configuration.

Method	Nodes	Sizes (Mvar)	Annual Costs (USD/year)	Red. (%)
Ben. Case	—	—	130,613.90	0.00
GNDO (Fixed)	$\left[\begin{array}{ccc}14& 25& 30\end{array}\right]$	$\left[\begin{array}{ccc}0.23083& 0.09996& 0.53905\end{array}\right]$	108,215.94	17.15
GNDO (Variable)	$\left[\begin{array}{ccc}14& 25& 30\end{array}\right]$	$\left[\begin{array}{ccc}0.24092& 0.10118& 0.69257\end{array}\right]$	106,550.69	18.42

compares fixed and variable reactive power injection by fixing the nodes reached by the GNDO approach in the exact MINLP model (1)–(8) using the BONMIN solver as a solution method.

The results in Table 6 show that an additional reduction of USD 1665.25 is obtained when variable reactive power injection is implemented with D-STATCOMs. That is to say, this is in comparison with the fixed injection scenario. This additional profit is reached when the total installed capabilities of the D-STATCOMs increase to 1034.67 kvar, which in turn implies an increase of 164.83 kvar with respect to the fixed injection scenario. With this additional reactive power, the D-STATCOMS have more flexibility to generate variable reactive power, i.e., they gain the ability to follow the demand behavior in order to inject only the needed reactive power for each hour as a function of the total grid requirements.

5.2. Meshed Configuration Results

Table 7. Optimal solutions reached by the comparison and proposed optimization methods with fixed reactive power injection in a meshed grid configuration.

Method	Nodes	Sizes (Mvar)	Annual Costs (USD/year)	Red. (%)
Ben. Case	—	—	86,914.74	0.00
COUENNE	$\left[\begin{array}{ccc}5& 6& 11\end{array}\right]$	$\left[\begin{array}{ccc}0.0000& 0.32577& 0.23872\end{array}\right]$	83,254.95	4.21
BONMIN	$\left[\begin{array}{ccc}15& 16& 17\end{array}\right]$	$\left[\begin{array}{ccc}0.11818& 0.00391& 0.34700\end{array}\right]$	81,171.76	6.61
SSA	$\left[\begin{array}{ccc}14& 30& 32\end{array}\right]$	$\left[\begin{array}{ccc}0.14620& 0.39440& 0.20230\end{array}\right]$	77,870.17	10.41
GNDO	$\left[\begin{array}{ccc}14& 30& 32\end{array}\right]$	$\left[\begin{array}{ccc}0.11554& 0.46482& 0.15147\end{array}\right]$	77,834.42	10.45

presents the comparison between the SSA approach, the BONMIN and COUENNE solvers, and the proposed GNDO approach that considers fixed reactive power injection throughout the day for the IEEE 33-node grid with a meshed grid topology.

The numerical results shown in Table 7 imply that:

• Once again, the COUENNE and BONMIN solvers got stuck in locally optimal solutions, which is attributed to the nonlinear non-convex nature of the exact MINLP model. The COUENNE solver only reached a reduction of $4.21\%$ with respect to the benchmark case, while the BONMIN solver yielded a better local solution, with a reduction of $6.61\%$ with respect to the benchmark.
• The proposed GNDO approach found a solution with an expected reduction in the annual operative costs of about 10.45%, i.e., 0.04% better than the solution reported by the SSA approach. However, the main characteristic of both solutions is that these located the D-STATCOMs in the same set of nodes (14, 30, and 32). However, their sizes differ, which explains the difference in the final objective function value.
• The total installed size of the D-STATCOMs with the GNDO is 731.83 kvar, whereas the SSA installed 742.90 kvar, which implies that with a better definition of the D-STATCOM sizes it is also possible to reach better final objective function values. Nevertheless, the main characteristic of this behavior is that based on its probability functions, the GNDO explored and exploited the solution space with better sensitivity than the SSA approach.

To confirm that the variable reactive power injections play an important role in reducing the expected annual operating costs of the network with respect to the fixed injection case, the results presented in

Table 8. Additional improvements reached when variable reactive power is used for daily compensation with D-STATCOMs for the meshed grid configuration.

Method	Nodes	Sizes (Mvar)	Annual Costs (USD/year)	Red. (%)
Ben. Case	—	—	86,914.74	0.00
GNDO (Fixed)	$\left[\begin{array}{ccc}14& 30& 32\end{array}\right]$	$\left[\begin{array}{ccc}0.11554& 0.46482& 0.15147\end{array}\right]$	77,834.42	10.45
GNDO (Variable)	$\left[\begin{array}{ccc}14& 30& 32\end{array}\right]$	$\left[\begin{array}{ccc}0.11578& 0.49915& 0.17417\end{array}\right]$	77.697.53	10.60

are employed.

The results in Table 8 show that the scenario with variable reactive power injection finds additional improvements regarding the final objective function value, i.e., approximately USD 136.89 per year of operation when compared with the fixed reactive power injection case. This result confirms that the main advantage of using D-STATCOMs is that variable reactive power is injected as a function of the network requirements, which justifies the distribution of static compensators instead of capacitor banks in order to support reactive power in distribution grids with variable demand behavior.

5.3. Additional Results

With the purpose of demonstrating the effectiveness and robustness of the GNDO approach in conjunction with the successive approximations power flow method, after 100 consecutive evaluations for the radial and meshed grid structures, the results in

Table 9. Statistical analysis for the radial and meshed configurations of the IEEE 33-bus grid.

Case	Min. (USD/year)	Max. (USD/year)	Mean (USD/year)	Std. Dev. (USD/year)	Time (s)
Radial	108,215.94	116,992.53	109,373.20	2053.04	42.70
Meshed	77,834.42	80,231.64	78,256.17	597.01	42.94

show the statistical analysis for both simulation cases.

The results in Table 9 show that:

• The difference between the minimum and maximum values reached in the radial configuration was about USD 8776.59, which implies that in the worst simulation case, the expected reductions with respect to the benchmark case would be 10.43%. When compared to the results shown in Table 5, this is better than both GAMS solutions. Even if the maximum result (worst result of the GNDO method) for the radial simulation case is a local solution, it has better characteristics than the local optimal solutions found with the MINLP BONMIN and COUENNE solvers in GAMS.
• As for the meshed configuration, the difference between the maximum and minimum objective function values was USD 2397.22, i.e., the worst solution reached by the GNDO approach reduces the expected annual operating costs by about 7.69%. This result confirms that for the meshed configuration case, the worst solution of the GNDO approach is also better than both solutions found with the GAMS software and with the BONMIN and COUENNE solvers.
• The mean values for the radial and meshed configurations are very close to the minimum values, i.e., most of the solutions are closer to each other in a closed ball with a radius equivalent to the standard deviation. These results confirm the effectiveness of the GNDO approach at solving the complex optimization problem involving the optimal location and size of D-STATCOMs in distribution grids via combinatorial optimization, with the main advantage that less than 43 s is required to reach a solution in both simulation cases.

6. Conclusions and Future Work

The optimal reactive power compensation problem via D-STATCOM location and sizing was studied in this research through the application of master–slave optimization methodology that solves its exact MINLP formulation. In the master stage, the GNDO approach was implemented using a discrete codification that decided the nodes where the D-STATCOMs were to be located in conjunction with their optimal sizes. In the slave stage, a power flow method (the successive approximations approach) was used to define the expected energy loss costs for a daily operation horizon. An improvement stage was implemented in the exact MINLP model that involved fixing the location of the D-STATCOMs provided by the master–slave optimizer in order to obtain an NLP model that could vary the reactive power injection on an hourly basis in order to obtain additional profits for the distribution company.

The numerical simulations run in the IEEE 33-bus system with radial and meshed configurations while including residential, industrial, and commercial users demonstrated that:

• In the radial configuration scenario, the proposed GNDO approach reached a lower minimum objective function value than the SSA approach. In addition, the improvement stage allowed for an additional reduction of USD 1665.25 regarding the fixed injection case, which confirms that the variable reactive power scenario via the daily control of the D-STATCOMs allows the distribution company to obtain additional profits in the total annual expected operating costs of the network.
• In the meshed configuration, the proposed GNDO approach found the same nodal location for the D-STATCOMs as the SSA methodology. However, the sizing was better, which allowed for a better final objective function value. In the variable reactive power scenario, an additional profit of USD 186.39 with respect to the fixed injection scenario was also reached, thus confirming that controlling the reactive power injection as a function of the system’s requirements is the better option to operate D-STATCOMs based on utilities.
• The expected annual operative gains for the meshed topology are considerably lower than those of the radial grid, i.e., USD 9080.32 for the meshed cased vs. USD 22,397.96 for the radial case. This is an expected result since in meshed configurations the energy losses are lower when compared to radial configurations, which is attributed to better flow distribution and voltage profiles.
• Regarding the processing times, the proposed GNDO approach found the solution for the MINLP model in less than one minute for both topologies, which confirms the effectiveness and robustness of the master–slave optimization methodology to deal with complex electrical engineering problems involving solution spaces with infinite dimensions and complex nonlinear, non-convex constraints.

In future work, the following studies can be conducted: (i) combining the D-STATCOMs with renewable generators in radial and meshed distribution grids by defining their sizes and locations in order to reduce the grid investments and the operating and maintaining costs; and (ii) a comparative study between different reactive power compensators such as D-STATCOMs, Thyristor-controlled series compensators (TCSCs), and Static VAR Compensators (SVCs), as well as their impact on the expected annual grid operating costs.