An Approximate Mixed-Integer Convex Model to Reduce Annual Operating Costs in Radial Distribution Networks Using STATCOMs

Abstract

The problem of optimal siting and sizing of distribution static compensators (STATCOMs) is addressed in this research from the point of view of exact mathematical optimization. The exact mixed-integer nonlinear programming model (MINLP) is decoupled into two convex optimization sub-problems, named the location problem and the sizing problem. The location problem is addressed by relaxing the exact MINLP model, assuming that all the voltages are equal to $1\angle 0^{\circ }$, which allows obtaining a mixed-integer quadratic programming model as a function of the active and reactive power flows. The solution of this model provides the best set of nodes to locate all the STATCOMs. When all the nodes are selected, it solves the optimal reactive power problem through a second-order cone programming relaxation of the exact optimal power flow problem; the solution of the SOCP model provides the optimal sizes of the STATCOMs. Finally, it refines the exact objective function value due to the intrinsic non-convexities associated with the costs of the STATCOMs that were relaxed through the application of Taylor’s series expansion in the location and sizing stages. The numerical results in the IEEE 33- and 69-bus systems demonstrate the effectiveness and robustness of the proposed optimization problem when compared with large-scale MINLP solvers in GAMS and the discrete-continuous version of the vortex search algorithm (DCVSA) recently reported in the current literature. With respect to the benchmark cases of the test feeders, the proposed approach reaches the best reductions with $14.17\%$ and $15.79\%$ in the annual operative costs, which improves the solutions of the DCVSA, which are $13.71\%$ and $15.30\%$, respectively.

1. Introduction

Electrical distribution networks represent the largest portion of power systems, which are entrusted with providing electrical energy from a transmission/sub-transmission substation to all end-users at medium- and low-voltage levels [1,2]. Owing to economic restrictions, investment distribution networks are typically constructed with radial configurations since this allows reducing inversion in the infrastructure, simplifying also the protective devices’ coordination [3]. However, the radial structure of these systems implies that they experience high power and energy losses when compared with transmission systems, being that the energy losses in these grids are between 6% and 18% of the total power input at the substation bus, while in transmissions networks, these are between 1.5% and 2.0% [4,5].

Owing to the high energy losses in distribution networks, utilities implement different methodologies to reduce them as much as possible. Some recurrent strategies are the optimal location of fixed-step capacitor banks [6,7,8]; optimal grid reconfiguration [9,10]; optimal siting and sizing of dispersed generators [11,12]; optimal siting and sizing of battery energy storage systems [13,14]; and reactive power compensation via distribution static compensators, i.e., STATCOMs [15,16,17]. In the case of fixed-step capacitor banks, these are attractive technologies due to their low cost and high reliability; however, these devices inject reactive power in discrete steps, which implies that, due to the daily load variability, the total grid power losses are not completely minimized. The strategies that use the grid reconfiguration to minimize power losses are efficient regarding the percentage of energy losses reduction; however, these require important modifications of the grid topology, including new line construction and redesign of the protective devices schemes, which is more expensive when compared with shunt reactive power compensation. The optimal integration of dispersed generators based on renewables is undoubtedly the most efficient option to reduce the total energy losses of the network, especially when active and reactive power capabilities of these devices are used [12]. However, from an economical point of view, the main application of the renewables is a reduction in greenhouse gas emissions and the reduction in the total energy purchasing costs in the substation buses. In the case of batteries, these devices help with the total energy losses reduction. However, their investment costs are bigger when compared with shunt reactive power compensation; the main application of the batteries is to extend the usability of the renewables in periods of time with high demand and low generation [18]. Finally, STATCOMs are excellent devices to reduce the amount of grid energy losses costs, with the main advantage that these can inject variable reactive power as a function of the grid load behavior, even if these are more expensive that the capacitor banks; these allow important improvements in the total annual grid operation costs. Moreover, these devices present useful life times of about 25 years, with minimum maintenance costs, which make these appropriate devices for improving the grid performance regarding energy losses and voltage profiles [19].

In the current literature regarding STATCOMs and their usage in distribution systems applications, some of the most relevant works are discussed below. Authors of [16] presented the application of the discrete-continuous version of the vortex-search algorithm (DCVSA), where the discrete part is entrusted with determining the nodes where the STATCOMs must be located and the continuous part defines their optimal sizes. Numerical results demonstrated that the DCVSA approach reaches the best optimal solution in the IEEE 33- and 69-bus systems when compared with the GAMS optimization package. Authors of [17] presented a complete revision of the strategies used to locate and size STATCOMs in distribution grids by making a classification in five categories, including analytical methods, artificial neural network-based approaches, metaheuristic methods, sensitivity approaches, and a combination of sensitivity approaches and metaheuristic methods. The main contribution of this research is the revision of the complete revision of the state of the art regarding distribution static compensators and their applications to improve the distribution system performance. In [20] was proposed a hybrid optimization approach based on the classical Chu and Beasley genetic algorithm in the master optimization stage that determines the nodes where the STATCOMs must be installed, while the slave stage uses a second-order cone programming reformulation of the multiperiod optimal power flow problem to determine their optimal sizes. Numerical results demonstrate the efficiency of the proposed approach when compared with the GAMS solvers in the IEEE 33- and 69-bus systems. Castiblanco-Pŕez et al. in [4] proposed the application of the discrete-continuous version of the Chu and Beasley genetic algorithm to locate and size STACOMs in distribution networks, with numerical results comparable with the DCVSA reported in [16]. The main contribution of the authors was the usage of differentiated load zones classified in residential, commercial, and industrial users, including the possibility of having radial and mesh distribution system topologies. Computational validations were made in the IEEE 33-bus systems, and numerical results demonstrated the possibility of having high-quality solutions when compared with the MINLP solvers available in the GAMS optimization package.

Even if the previous optimization methodologies are efficient to solve the problem addressed in this research, these have two main problems, which are (i) the usage of combinatorial methods to locate and size the STATCOMs since their random exploration and exploitative nature through the solution space makes it impossible to ensure the global optimum finding, and (ii) in the case of the exact MINLP approaches, the main complication is associated with the non-convexity of the solution space in their continuous component that highly increases the possibility of staying stuck in local optimal solutions. To deal with these problems, this research proposes a new methodology based on the hybridization of two convex-based models. The proposed methodology is completely described below.

Unlike aforementioned works regarding the optimal siting and dimensioning of STATCOMs in distribution networks, in this research, we propose a new optimization methodology based on mixed-integer convex optimization that decouples the location problem from the sizing problem. To solve the location problem, a mixed-integer quadratic formulation is employed considering that the voltage profiles are equal to $1.0\angle 0^{\circ }$. This becomes the power balance equations in a set of linear equations that becomes the location problem into a mixed-integer approximation that ensures the optimality of the solution. The solution of the location problem allows fixing the values of the binary variables in the sizing problem. The assignation of the binary variables allow solving the sizing problem of the STATCOMs via the second-order cone programming reformulation proposed in [20] with the main advantage that the global optimum finding is ensured via convex optimization. Numerical results in the IEEE 33- and 69-bus systems demonstrate the effectiveness and robustness of the proposed optimization approach when compared with the DCVSA and the MINLP solvers in GAMS.

The main contributions of this research are the following:

• The presentation of a new optimization methodology to solve separately the problems of optimal location and sizing of STATCOMs in distribution networks to reduce the annual operative costs of the network, where the location problem is solved through a mixed-integer quadratic formulation, and the sizing stage is addressed with a second-order cone programming equivalent.
• The validation of the proposed methodology in two classical test feeders composed of 33 and 69 nodes with better results than the best current approach reported in the current literature, i.e., the discrete-continuous version of the vortex search algorithm.

It is worth mentioning that in the scope of this research, distribution networks with meshed structures are not considered since the second-order cone programming reformulation of the optimal power flow problem proposed in [21] only deals with purely-radial distribution grids; however, this can be taken as a opportunity of research in future works.

The proposed decoupled optimization approaches to solve the problem of the location and sizing of STATCOMs in distribution networks based on mixed-integer convex and conic programming was not previously proposed in the current literature. This was identified as a gap in the current literature that this research work is aiming to fulfill.

This document is structured as follows: Section 2 presents the general optimization problem regarding the optimal placement and sizing of STATCOMs in distribution grids with a radial structure. Section 3 presents the main aspects of the proposed solution methodology based on decoupling the location from the sizing problem via mixed-integer quadratic and second-order cone programming. Section 4 shows the main characteristics of the IEEE 33- and 69-bus systems, including the parametrization of the costs function of the STATCOMs. Section 5 presents the main numerical results of the proposed solution methodology and their comparison with recent literature reports. Section 6 lists the main conclusions obtained from this research and defines some possible future works.

2. Optimization Problem

The problem of the optimal siting and sizing of STATCOMs in distribution networks corresponds to a mixed-integer nonlinear programming (MINLP) problem that involves binary variables associated with the location or not of a STATCOM in a particular node, and continuous variables regarding power flow equations, voltages and currents, among others. The complete MINLP model is presented below.

2.1. Objective Function

The objective function for the problem of the optimal annual cost reduction in distribution networks considering STATCOMs corresponds to the minimization of the energy losses cost in conjunction with the total investment cost. The objective function can be formulated as follows:

$$\begin{array}{c}\hfill min{A}_{\mathrm{cost}}=z_1+z_2,\end{array}$$

(1)

$$\begin{array}{c}\hfill z_1=C_{\mathrm{kWh}}T\sum _{ij\in \mathcal{L}}\sum _{h\in \mathcal{H}}{R}_{ij}{i}_{ij,h}^2\Delta h,\end{array}$$

(2)

$$\begin{array}{c}\hfill z_2=T\left(\frac{k_1}{k_2}\right)\sum _{j\in \mathcal{N}}\left(\alpha y_j^2+\beta y_j+\gamma \right)y_j,\end{array}$$

(3)

where $A_{\mathrm{cost}}$ is the objective function value regarding the summation of the annual operating and investment costs; $z_1$ is the value of the annual operative costs of the grid associated with the energy losses; $z_2$ is a cubic function regarding the investment costs in STATCOMs; $C_{\mathrm{kWh}}$ is the average cost of the energy losses; T is the number of days in an ordinary year (i.e., 365 days); $R_{ij}$ is the resistance parameter associated with the line that connects nodes i and j; $i_{ij}$ is the current magnitude that flows through this line; $\Delta h$ is the length of the period of time where electrical variables are constant (typically 0.25 h, 0.50 h, or 1 h); $k_1$ and $k_2$ are parameters associated with the annualization of investment costs of the STATCOMs; and $\alpha$, $\beta$, and $\gamma$ are the cubic, quadratic, and linear costs coefficients associated with the size $y_j$ assigned to the STATCOM connected at node j. Note that $\mathcal{L}$, $\mathcal{H}$, and $\mathcal{N}$, are the sets that contain all the lines, periods of time, and nodes of the networks, respectively.

Remark 1.

The component $z_1$ of the objective function is a convex quadratic function; however, the component $z_2$ is a cubic non-convex function that makes that the total annual operative costs take a non-convex structure.

2.2. Set of Constraints

The set of constraints in the problem of the optimal placement and sizing of STATCOMs in distribution networks includes the active and reactive power balance equations, the voltage drops at each line, and the maximum and minimum voltage bounds and device capabilities among others. The complete list of constraints is listed below.

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

(4)

$$\begin{array}{c}\hfill q_{ij,h}-X_{ij}{i}_{ij,h}^2-\sum _{k:\left(jk\right)\in \mathcal{L}}{q}_{jk,h}+q_{j,h}^{st}=Q_{j,h}^d,\left\{\forall j\in \mathcal{N},j\ne \mathrm{slack},\forall h\in \mathcal{H}\right\},\end{array}$$

(5)

$$\begin{array}{c}\hfill v_{j,h}^2=v_{i,h}^2-2\left(R_{ij}{p}_{ij,h}+X_{ij}{q}_{ij,h}\right)+\left(R_{ij}^2+X_{ij}^2\right)i_{ij,h}^2\left\{\forall ij\in \mathcal{E},\forall h\in \mathcal{H}\right\},\end{array}$$

(6)

$$\begin{array}{c}\hfill p_{ij,h}^2+q_{ij,h}^2=v_{i,h}^2{i}_{ij,h}^2,\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\}\end{array}$$

(7)

$$\begin{array}{c}\hfill x_j{q}_j^{st,min}\le y_j\le x_j{q}_j^{st,max}\left\{\forall j\in \mathcal{N}\right\},\end{array}$$

(8)

$$\begin{array}{c}\hfill -y_j\le q_{j,h}^{st}\le y_j\left\{\forall j\in \mathcal{N},\forall h\in \mathcal{H}\right\},\end{array}$$

(9)

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

(10)

$$\begin{array}{c}\hfill -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}$$

(11)

$$\begin{array}{c}\hfill \sum _{i\in \mathcal{N}}{x}_j\le N_{st}^{max},\end{array}$$

(12)

$$\begin{array}{c}\hfill x_j\in \left\{0,1\right\},\left\{\forall i\in \mathcal{N}\right\},\end{array}$$

(13)

where $p_{ij,h}$$\left(p_{jk,h}\right)$ and $q_{ij,h}$$\left(p_{jk,h}\right)$ are the active and reactive power flows through the line $ij$$\left(jk\right)$ at the period of time h; $P_{j,h}^d$ and $Q_{j,h}^d$ represent the active and reactive power consumptions at node j in the period of time h; $q_{j,h}^{st}$ is the reactive power injection in the STATCOM connected at node j at each period of time; $X_{ij}$ is the reactance parameter of the line between nodes i and j; $v_{i,h}$ and $v_{j,h}$ are the voltage magnitudes at nodes i and j in the period of time h; $q_j^{st,min}$ and $q_j^{st,max}$ are the minimum and maximum reactive power injection bounds allowed for a STATCOM connected at node j; $x_j$ is the binary variable that decides if a STATCOM is connected ($x_j=1$) or not ($x_j=0$) in the node j; $v_j^{min}$ and $v_j^{max}$ are the minimum and maximum voltage regulation bounds permitted by regulatory policies in all the nodes of the network; $i_{ij}^{max}$ is the thermal bound (maximum current) associated with the conductor that connects nodes i and j; and $N_{st}^{max}$ is the maximum number of STATCOMs available for installation along the distribution grid.

2.3. Model Interpretation

The complete interpretation of the mathematical model defined from (1) to (13) is the following: Equation (1) is the objective function value associated with the summation of the operative costs regarding the annual energy losses, i.e., Equation (2), and the annual investment costs regarding the installation of the STATCOMs along the grid, i.e., Equation (3). Equality constraints (4) and (5) represents the active and reactive power equilibrium equations in all nodes of the network at each period of time, respectively; Equation (6) represents the voltage drop at each line as a function of the current and active and reactive power flows; Equation (7) defines the application of Tellegen’s theorem at each sending power flow through the lines of the network at each period of time; and box-type constraints (8) and (9) define the possibility of installation of a STATCOM at node j and its maximum and minimum reactive power injection capabilities, respectively. Box-type constraints (10) and (11) define the grid voltage regulation limits and the maximum capabilities of the lines to transport current (i.e., thermal capabilities of the lines), respectively. Inequality constraint (12) limits the maximum number of STATCOMs that can be installed in the whole distribution network. Constraint (13) shows the binary nature of the decision variable.

Remark 2.

The main characteristic of the set of constraints (4)–(13) is the non-convexity of the solution space, specially caused by power equilibrium constraints, the voltage drops and the application of Tellegen’s theorem to define the hyperbolic relation between voltages, currents and powers at each line, respectively.

Owing to the complication of the optimization model, (1)–(13), that defines the optimal placement and sizing of STATCOMs in distribution grids with radial structure, this research proposes a mixed-integer convex approximation to solve the optimization by decoupling the optimal location problem from the optimal dimensioning problem. The proposed methodology is completely presented in the next section.

3. Solution Methodology

To deal with the problem of the optimal siting and sizing of STATCOMs in distribution networks, this research proposes a methodology based on the decoupling of the problem in the siting problem and sizing problem, respectively. In the problem of optimal sizing, we propose a mixed-integer quadratic that decides the best set of nodes for locating STATCOMs, and the sizing problem is solved in the MINLP model (1)–(13) by reducing it to an NLP equivalent.

3.1. Solution of the Location Problem

The location problem of the STATCOMs in distribution networks corresponds to the part of the complete optimization model associated with the definition of the binary variables, i.e., $x_j$. To deal with this part of the problem, we employ the mixed-integer quadratic approximation proposed in [22] to solve the optimal reconfiguration problem in pure-radial distribution networks. The author of [22] proposes that the power flow equations in the optimization model (1)–(13) can be reduced by assuming that all the voltages are near the unity in per-unit representation; moreover, we can assume that the resistance and reactance power losses in the lines are negligible with respect to the magnitudes of the active and reactive power flows. With this simplification, the set of constraints (4)–(13) can be reduced as follows:

$$\begin{array}{c}\hfill p_{ij,h}-\sum _{k:\left(jk\right)\in \mathcal{L}}{p}_{jk,h}=P_{j,h}^d,\left\{\forall j\in \mathcal{N},j\ne \mathrm{slack},\forall h\in \mathcal{H}\right\},\end{array}$$

(14)

$$\begin{array}{c}\hfill q_{ij,h}-\sum _{k:\left(jk\right)\in \mathcal{L}}{q}_{jk,h}+q_{j,h}^{st}=Q_{j,h}^d,\left\{\forall j\in \mathcal{N},j\ne \mathrm{slack},\forall h\in \mathcal{H}\right\},\end{array}$$

(15)

$$\begin{array}{c}\hfill x_j{q}_j^{st,min}\le y_j\le x_j{q}_j^{st,max}\left\{\forall j\in \mathcal{N}\right\},\end{array}$$

(16)

$$\begin{array}{c}\hfill -y_j\le q_{j,h}^{st}\le y_j\left\{\forall j\in \mathcal{N},\forall h\in \mathcal{H}\right\},\end{array}$$

(17)

$$\begin{array}{c}\hfill \sum _{i\in \mathcal{N}}{x}_j\le N_{st}^{max},\end{array}$$

(18)

$$\begin{array}{c}\hfill x_j\in \left\{0,1\right\},\left\{\forall i\in \mathcal{N}\right\},\end{array}$$

(19)

The main advantage of the set of constraints (14)–(19) is that it is a mixed-integer linear group of constraints that can be solvable with a modified version of the branch and bound method. Moreover, component $z_1$ of the objective function can be replaced with the help of (7) as follows (this replacement can be made if we assume that all the voltages are equal to $1.0\angle 0^{\circ }$):

$$\begin{array}{c}\hfill z_1^{\mathrm{approx}}=C_{\mathrm{kWh}}T\sum _{ij\in \mathcal{L}}\sum _{h\in \mathcal{H}}{R}_{ij}\left(p_{ij,h}^2+q_{ij,h}^2\right)\Delta h.\end{array}$$

(20)

In the case of the annual investment costs in STATCOMs, this can be approximated with the Taylor’s series expansion in the point $y_j=0$ since for values between $0$ Mvar and $10$ Mvar, their tendency is basically linear. With this approximation, the component $z_2$ of the objective function (1) can be approximated as presented below:

$$\begin{array}{c}\hfill z_2^{\mathrm{approx}}=T\gamma \left(\frac{k_1}{k_2}\right)\sum _{j\in \mathcal{N}}{y}_j.\end{array}$$

(21)

Remark 3.

Note that the relaxed set of constraints (14)–(19) with the algebraic sum of components $z_1^{\mathrm{approx}}$ and $z_2^{\mathrm{approx}}$ in (20) and (21) produce a mixed-integer quadratic convex model, which can be solved by ensuring the global optimum finding with the appropriate combination of the branch and bound method with interior point methods [23].

3.2. Solution of the Sizing Problem

Once the binary variables are obtained from the solution of the location problem presented in the previous section, the set of variables $x_j$ is fixed in the optimization model (1)–(13), which changes this model from an MINLP problem into a NLP problem. However, due to the non-convexity of the solution space, it is not possible to ensure the global optimum finding; for this reason, as recommended in [20], the equivalent second-order cone programming model can be obtained for it if we define two auxiliary variables, $u_{j,h}=v_{j,h}^2$ and $l_{ij,h}=i_{ij,h}^2$. The second-order cone equivalent set of constraints for the problem of the optimal sizing of STATCOMs in pure-radial distribution networks takes the following form:

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

(22)

$$\begin{array}{c}\hfill q_{ij,h}-X_{ij}{i}_{ij,h}-\sum _{k:\left(jk\right)\in \mathcal{L}}{q}_{jk,h}+q_{j,h}^{st}=Q_{j,h}^d,\left\{\forall j\in \mathcal{N},j\ne \mathrm{slack},\forall h\in \mathcal{H}\right\},\end{array}$$

(23)

$$\begin{array}{c}\hfill u_{j,h}=u_{i,h}-2\left(R_{ij}{p}_{ij,h}+X_{ij}{q}_{ij,h}\right)+\left(R_{ij}^2+X_{ij}^2\right)i_{ij,h}\left\{\forall ij\in \mathcal{E},\forall h\in \mathcal{H}\right\},\end{array}$$

(24)

$$\begin{array}{c}\hfill p_{ij,h}^2+q_{ij,h}^2+\frac{1}{4}{\left(u_{i,h}-l_{ij,h}\right)}^2\le \frac{1}{4}{\left(u_{i,h}+l_{ij,h}\right)}^2,\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\}\end{array}$$

(25)

$$\begin{array}{c}\hfill x_j{q}_j^{st,min}\le y_j\le x_j{q}_j^{st,max}\left\{\forall j\in \mathcal{N}\right\},\end{array}$$

(26)

$$\begin{array}{c}\hfill -y_j\le q_{j,h}^{st}\le y_j\left\{\forall j\in \mathcal{N},\forall h\in \mathcal{H}\right\},\end{array}$$

(27)

$$\begin{array}{c}\hfill {\left(v_j^{min}\right)}^2\le u_{j,h}\le {\left(v_j^{min}\right)}^2\left\{\forall j\in \mathcal{N},\forall h\in \mathcal{H}\right\},\end{array}$$

(28)

$$\begin{array}{c}\hfill l_{ij,h}\le {\left(i_{ij}^{max}\right)}^2\left\{\forall ij\in \mathcal{L},\forall h\in \mathcal{H}\right\},\end{array}$$

(29)

Remark 4.

Note that the constraint (25) is obtained if the product of two continuous variables $u_{i,h}{l}_{ij,h}$ is replaced by its equivalent hyperbolic formulation, i.e.,

$$\begin{array}{c}\hfill 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\end{array}$$

In the case of the objective function, to ensure that it is also convex, we approximate $z_2$ as presented in (21), which allows obtaining the following approximate annual cost of the network:

$$\begin{array}{c}\hfill min{A}_{\mathrm{cost}}=z_1+z_2^{\mathrm{approx}},\end{array}$$

(30)

$$\begin{array}{c}\hfill z_1=C_{\mathrm{kWh}}T\sum _{ij\in \mathcal{L}}\sum _{h\in \mathcal{H}}{R}_{ij}{l}_{ij,h}\Delta h,\end{array}$$

(31)

$$\begin{array}{c}\hfill z_2^{\mathrm{approx}}=T\gamma \left(\frac{k_1}{k_2}\right)\sum _{j\in \mathcal{N}}{y}_j,\end{array}$$

(32)

Finally, the solution of the optimization model (22)–(32) ensures the global optimum finding of the sizes of the STATCOMs due to the convexity of the solution space [24].

3.3. Methodology Summary

The proposed optimization methodology to solve the problem of the optimal placement and sizing of STATCOMs in distribution networks is summarized in the flowchart diagram depicted in Figure 1.

Note that the proposed optimization methodology depicted in Figure 1 is independent of the optimization tool and it can be implemented in multiple software that deal with mixed-integer convex formulations, with the main advantage that its solution will always be the same due to the convexity of the solution space in the case that the binary variable input in the second optimization model is the same.

4. Test Feeders

The computational validation of the proposed decoupled methodology to determine the optimal size and location of STATCOMs in distribution networks is carried out in two classical test feeders composed of 33 and 69 buses. These test feeders are depicted in Figure 2.

The main characteristic of these test feeders is that both operate at the substation bus with an output voltage of 12.66 kV. The electrical parameters for these grids are reported in

Table 1. Branch and load parameters of the IEEE 33-bus system.

Node i	Node j	${\mathit{R}}_{\mathit{ij}}$ ($\Omega$)	${\mathit{X}}_{\mathit{ij}}$ ($\Omega$)	${\mathit{P}}_{\mathit{j}}$ (kW)	${\mathit{Q}}_{\mathit{j}}$ (kvar)	Node i	Node j	${\mathit{R}}_{\mathit{ij}}$ ($\Omega$)	${\mathit{X}}_{\mathit{ij}}$ ($\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

and

Table 2. Branch and load parameters of the IEEE 69-bus system.

Node i	Node j	${\mathit{R}}_{\mathit{ij}}$ ($\Omega$)	${\mathit{X}}_{\mathit{ij}}$ ($\Omega$)	${\mathit{P}}_{\mathit{j}}$ (kW)	${\mathit{Q}}_{\mathit{j}}$ (kvar)	Node i	Node j	${\mathit{R}}_{\mathit{ij}}$ ($\Omega$)	${\mathit{X}}_{\mathit{ij}}$ ($\Omega$)	${\mathit{P}}_{\mathit{j}}$ (kW)	${\mathit{Q}}_{\mathit{j}}$ (kvar)
1	2	0.0005	0.0012	0	0	3	36	0.0044	0.0108	26	18.55
2	3	0.0005	0.0012	0	0	36	37	0.0640	0.1565	26	18.55
3	4	0.0015	0.0036	0	0	37	38	0.1053	0.1230	0	0
4	5	0.0251	0.0294	0	0	38	39	0.0304	0.0355	24	17
5	6	0.3660	0.1864	2.6	2.2	39	40	0.0018	0.0021	24	17
6	7	0.3810	0.1941	40.4	30	40	41	0.7283	0.8509	1.2	1
7	8	0.0922	0.0470	75	54	41	42	0.3100	0.3623	0	0
8	9	0.0493	0.0251	30	22	42	43	0.0410	0.0475	6	4.3
9	10	0.8190	0.2707	28	19	43	44	0.0092	0.0116	0	0
10	11	0.1872	0.0619	145	104	44	45	0.1089	0.1373	39.22	26.3
11	12	0.7114	0.2351	145	104	45	46	0.0009	0.0012	39.22	26.3
12	13	1.0300	0.3400	8	5	4	47	0.0034	0.0084	0	0
13	14	1.0440	0.3450	8	5.5	47	48	0.0851	0.2083	79	56.4
14	15	1.0580	0.3496	0	0	48	49	0.2898	0.7091	384.7	274.5
15	16	0.1966	0.0650	45.5	30	49	50	0.0822	0.2011	384.7	274.5
16	17	0.3744	0.1238	60	35	8	51	0.0928	0.0473	40.5	28.3
17	18	0.0047	0.0016	60	35	51	52	0.3319	0.1114	3.6	2.7
18	19	0.3276	0.1083	0	0	9	53	0.1740	0.0886	4.35	3.5
19	20	0.2106	0.0690	1	0.6	53	54	0.2030	0.1034	26.4	19
20	21	0.3416	0.1129	114	81	54	55	0.2842	0.1447	24	17.2
21	22	0.0140	0.0046	5	3.5	55	56	0.2813	0.1433	0	0
22	23	0.1591	0.0526	0	0	56	57	1.5900	0.5337	0	0
23	24	0.3460	0.1145	28	20	57	58	0.7837	0.2630	0	0
24	25	0.7488	0.2475	0	0	58	59	0.3042	0.1006	100	72
25	26	0.3089	0.1021	14	10	59	60	0.3861	0.1172	0	0
26	27	0.1732	0.0572	14	10	60	61	0.5075	0.2585	1244	888
3	28	0.0044	0.0108	26	18.6	61	62	0.0974	0.0496	32	23
28	29	0.0640	0.1565	26	18.6	62	63	0.1450	0.0738	0	0
29	30	0.3978	0.1315	0	0	63	64	0.7105	0.3619	227	162
30	31	0.0702	0.0232	0	0	64	65	1.0410	0.5302	59	42
31	32	0.3510	0.1160	0	0	11	66	0.2012	0.0611	18	13
32	33	0.8390	0.2816	14	10	66	67	0.0047	0.0014	18	13
33	34	1.7080	0.5646	19.5	14	12	68	0.7394	0.2444	28	20
34	35	1.4740	0.4873	6	4	68	69	0.0047	0.0016	28	20

, respectively.

To emulate the daily behavior of the active and reactive power consumption of the network,

Table 3. Daily active and reactive demand profiles.

Period	Act. (pu)	React. (pu)	Period	Act. (pu)	React. (pu)
1	0.1700	0.1477	25	0.4700	0.3382
2	0.1400	0.1119	26	0.4700	0.3614
3	0.1100	0.0982	27	0.4500	0.3877
4	0.1100	0.0833	28	0.4200	0.3434
5	0.1100	0.0739	29	0.4300	0.3771
6	0.1000	0.0827	30	0.4500	0.4269
7	0.0900	0.0831	31	0.4500	0.4224
8	0.0900	0.0637	32	0.4500	0.3647
9	0.0900	0.0702	33	0.4500	0.4226
10	0.1000	0.0875	34	0.4500	0.3081
11	0.1100	0.0728	35	0.4500	0.2994
12	0.1300	0.1214	36	0.4500	0.3336
13	0.1400	0.1231	37	0.4300	0.3543
14	0.1700	0.1390	38	0.4200	0.3399
15	0.2000	0.1410	39	0.4600	0.4234
16	0.2500	0.1998	40	0.5000	0.4061
17	0.3100	0.2497	41	0.4900	0.3820
18	0.3400	0.3224	42	0.4700	0.3820
19	0.3600	0.3263	43	0.4500	0.3887
20	0.3900	0.3661	44	0.4200	0.2751
21	0.4200	0.3585	45	0.3800	0.3383
22	0.4300	0.3316	46	0.3400	0.2355
23	0.4500	0.4187	47	0.2900	0.2301
24	0.4600	0.3652	48	0.2500	0.1818

reports the daily variation of these consumptions in steps of half hours. Note that these values must be scaled by 2 during all simulations [16].

The assessment of the component of the objective function $z_2$ can be made with the parameters reported in

Table 4. Parameters to evaluate the component $z_2$ of the objective function.

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

. Some of these values were adapted from [25] and complemented with the information provided in [16].

To evaluate the objective function defined in (1), the parameters reported in Table 4 are considered. Some of these parameters were taken from [16,25]. To avoid mistakes during the evaluation of this objective function, it is important to mention that the variable $y_j$ must be defined in Mvar.

5. Computational Validation

The solution of the MINLP model defined from (1)–(13) and the decoupled models (14)–(21) and (22)–(32) were implemented in GAMS software with the solvers BONMIN and COUENNE and with the CVX programming tool in MATLAB with the help of the Gurobi solver. In the case of the MATLAB implementations, we used $2021b$ on a PC with an AMD Ryzen 7 3700 $2.3$-GHz processor and 16.0 GB RAM, running on a 64-bit version of Microsoft Windows 10 Single language. Moreover, the numerical results are also compared with the discrete-continuous version of the vortex search algorithm (DCVSA) proposed in [16].

5.1. Results in the IEEE 33-Bus System

Table 5. Numerical results in the IEEE 33-bus system for the proposed and comparative methods.

Method	Location (node)	Size (kvar)	${\mathit{A}}_{\mathbf{cost}}$ (US$/year)
Benchmark case	—	—	112,740.90
BONMIN [26]	$\left\{18,19,33\right\}$	$\left\{171.15,0.00,472.69\right\}$	103,945.49
COUENNE [26]	$\left\{6,7,12\right\}$	$\left\{296.47,165.64,248.4\right\}$	98,936.36
DCVSA [16]	$\left\{14,30,32\right\}$	$\left\{159.90,359.10,107.20\right\}$	97,284.49
Proposed	$\left\{14,30,32\right\}$	$\left\{196.22,439.80,141.20\right\}$	96,767.69

presents the numerical comparison between the GAMS solvers, the DCVSA and the proposed decoupled methodology.

Numerical results in Table 5 show that: (i) the proposed decoupled methodology identifies nodes 14, 30, and 32 to locate all the STATCOMs, which coincides with the solution reported by the DCVSA. However, the objective function obtained by our proposed approach is USD 516.80 more economic than the DCVSA, which is attributable to the optimal sizes reached by the proposed conic model in the sizing stage, being that this is greater than the sizes provided by the DCVSA in [16]; (ii) the exact approaches in the GAMS software with the BONMIM and COUENNE solvers are stuck in local optimal solutions, which is attributable to the non-convexity of the exact model (1)–(13); and (iii) with respect to the benchmark case, the proposed approach reaches a reduction of $14.17\%$ followed by the DCVSA with $13.71\%$, and the COUENNE solution with $12.24\%$, respectively.

To illustrate the daily dynamic performance of the proposed approach regarding the reactive power injection in the selected nodes of the network, in Figure 3 are reported the reactive power injections in all the periods of time for all the STATCOMs. In general, it is observed that between the periods of time 1 and 17, all the STATCOMs control their reactive power injection below the nominal capacity, as during these time periods, the active and reactive power demand of the grid is in a general low; however, after period 18, all the STATCOMs inject their rate capacities since in these time periods, the grid demand increases considerably. Finally, at the end of the operation period, the reactive power injections in the STATCOMs start to reduce, owing to the reduction in the apparent power consumption on the whole grid.

The most important fact that Figure 3 tries to show is that the reactive power compensation in distribution networks is indeed a dynamic compensation problem since the grid requirements vary as a function of the total active and reactive demand during the day. In this sense, depending on the active and reactive demand curves employed, the optimal location and sizes of the STACOMS can vary significantly as shown in [4], where residential, industrial, and commercial users were considered along with the distribution feeders.

5.2. Results in the IEEE 69-Bus System

Table 6. Numerical results in the IEEE 69-bus system for the proposed and comparative methods.

Method	Location (node)	Size (kvar)	${\mathit{A}}_{\mathbf{cost}}$ (US$/year)
Benchmark case	—	—	119,715.63
DCVSA [16]	$\left\{21,61,64\right\}$	$\left\{83.87,460.06,113.88\right\}$	101,399.89
Proposed	$\left\{22,61,64\right\}$	$\left\{92.80,574.49,150.55\right\}$	100,806.91

shows the comparison between the DCVSA and the proposed decoupled methodology. Note that solutions with the BONMIN and the COUENNE solvers in the GAMS software are not provided since these failed in the case of the IEEE 69-bus system due to the complexity of the solution space and the increase in the number of decision variables.

Numerical results in Table 6 show that (i) the DCVSA is stuck in a local optimum since it is obtained as the optimal location of the set of nodes 21, 61 and 64, while the proposed approach identifies node 22 instead of node 21; (ii) with respect to the benchmark case, the proposed approach could reach an annual cost reduction of $15.79\%$, while the DCVSA obtains a reduction of $15.30\%$. These values imply an additional gain of USD 592.98 per year of operation; and (iii) the optimal sizes of the STATCOMs reached by our proposal are bigger than the sizes provided by the DCVSA. However, this situation happens due to the fact that the DCVSA works in the first stage with constant power injections during all periods of time to define the STATCOMs’ sizes; then, an optimal multiperiod power flow is applied to refine their operation, while our proposal works directly under a multiperiod programming scenario that allows identifying better solutions as reported in Table 6.

On the other hand, Figure 4 shows the dynamic reactive power injection in all the STATCOMs installed for the IEEE 69-bus system.

Results in Figure 4 show in general that the STATCOMs inject dynamic reactive power as a function of the grid requirements. Note that in the periods between 1 and 18, these STATCOMs work with reactive power injections below their nominal capacities; however, when the demand increases considerably, these STATCOMs inject their maximum rates (see periods after 18 h).

It is worth mentioning that the behavior in the IEEE 69-bus system is very similar to the IEEE 33-bus case; which is explainable by the fact that for both systems, the same daily active and reactive power load profiles reported in Table 3 are used.

Note that Figure 3 and Figure 4 confirm that the problem of the optimal reactive power compensation in distribution grids is indeed a dynamical optimization problem, where the final solution is highly influenced by the characterization of the demand curve of the grid in terminals of the substation. This shows the importance of developing a previous study to the application of the proposed decoupled optimization approach that corresponds to the determination of the average load profile of the distribution. It is noteworthy that this curve is the main external input for the application of any optimization methodology, even if these are based on metaheuristics or exact approaches.

6. Conclusions and Future Works

The problem of the optimal siting and sizing STATCOMs in radial distribution grids to reduce the annual grid operative costs was addressed in this research from the point of view of mixed-integer convex optimization. The exact MINLP model was decoupled into two subproblems, named the location problem and the sizing problem. The former problem was formulated with a mixed-integer quadratic programming model that allows determining the optimal location of the STATCOMs. These locations are provided in the latter stage (sizing problem), where a second-order cone programming model was used to solve the multiperiod optimal power flow problem. To obtain these equivalent models, the component of the objective function regarding the annualized costs of the STATCOMs was linearized through Taylor’s series expansion. Numerical results demonstrated that the proposed decoupled approach improves the solutions reached by the DCVSA and the BONMIN and COUENNE solvers available in the GAMS optimization package for both test feeders under study.

In the case of the IEEE 33-bus system, the general improvement with respect to the base case was $14.17\%$, which implied a reduction of about USD 15,973.21 per year of operation; for the IEEE 69-bus system, this reduction was $15.79\%$, i.e., USD 18,908.72 per year of operation. These results improved the best numerical solution reported in the literature with the DCVSA of about USD 516.80 for the IEEE 33-bus system and USD 592.98 in the case of the IEEE 69-bus system, which confirmed the effectiveness and robustness of the proposed approach regarding literature reports.

As future works, it will be possible to develop the following research studies: (i) to extend the proposed decoupled optimization methodology to integrate renewable sources (i.e., solar photovoltaic or wind sources) to minimize the total grid generation costs; (ii) to extend the proposed formulation to distribution grids with meshed topologies; and (iii) to apply the optimization paradigm known as poly-optimization to solve problems regarding the optimal placement and sizing of shunt devices in distribution networks.