Simultaneous Minimization of Energy Losses and Greenhouse Gas Emissions in AC Distribution Networks Using BESS

Abstract

The problem of the optimal operation of battery energy storage systems (BESSs) in AC grids is addressed in this paper from the point of view of multi-objective optimization. A nonlinear programming (NLP) model is presented to minimize the total emissions of contaminant gasses to the atmosphere and costs of daily energy losses simultaneously, considering the AC grid complete model. The BESSs are modeled with their linear relation between the state-of-charge and the active power injection/absorption. The Pareto front for the multi-objective optimization NLP model is reached through the general algebraic modeling system, i.e., GAMS, implementing the pondered optimization approach using weighting factors for each objective function. Numerical results in the IEEE 33-bus and IEEE 69-node test feeders demonstrate the multi-objective nature of this optimization problem and the multiple possibilities that allow the grid operators to carry out an efficient operation of their distribution networks when BESS and renewable energy resources are introduced.

1. Introduction

The problem of the optimal operation of battery energy storage systems (BESSs) in electrical distribution has attracted much interest in the last two decades due to the advances in energy storage and power electronic converter technologies [1,2,3]. The main advantages of including these batteries are as follows: (i) possible greenhouse emissions minimization in power systems fed by fossil generation [4,5]; (ii) improvement of the voltage profiles and the reduction of the grid power losses [6]; (iii) possible reduction of power oscillations caused by uncertainties in the renewable energy resources [7,8,9]; (iv) dynamic active and reactive power compensation using power electronic converters in electric vehicles [10]; and (v) real-time demand management [11]. The main question that must be solved when batteries are included in AC grids is as follows: How must these batteries operate to improve the electrical grid performance? Note that the solution of this question involves the proposition and solution of an optimal multi-period power flow model [12], which is known in the literature as the economic dispatch model when the objective function is concerned with operative cost minimization [2].

Nowadays, BESSs are real options for improving the electrical performance of distribution networks in medium- and low-voltage levels [1,13], since these devices can facilitate peak energy consumption fulfillment, enhance the profits from the integration of renewables and distributed energy sources, help power quality management, and generally reduce distribution network expansion costs [14]. The authors of [15] presented an experimental case with regards of the designing a BESS for urban industrial application in order to improve local power quality. They considered two cases of study, i.e., reactive power compensation and rapid variations of the electricity variations caused by a metro train load and a production of a solar power plant, and they mentioned that these studies were technically successful. In addition, an evaluation of the BESS assessment in a real distribution system in Riyadh, Saudi Arabia, was studied by Almehizia et al. [16]. Numerical results quantify the economic return that BESSs can offer for the electrical distribution grids, thus showing that improving the decision-making of the electric utility can deal with increasing load demand and power quality issues.

To analyze BESS in electrical networks with AC or DC operational technologies, one of the most accepted and currently used models corresponds to the linear representation between the state-of-charge and the power injection/absorption [17,18]. However, even if the battery model can be represented using linear equations, nonlinearities appear due to the nonlinear non-convex behaviors in AC power grids, which are caused by the power balance equations [1,19]. To address the nonlinear economic dispatch problem involving BESS and renewable generation in AC grids, various optimization models with different grades of approximation can be used. The first model corresponds to the adaption of the classical economic dispatch model for thermal power systems that uses a unique nodal representation of the system [18]; the second model is the DC approximation of the power flow equations that works with the angular difference between two adjacent nodes by assuming that the voltage is equal to 1.0 pu [20]; and the third model corresponds to the complete AC power flow model that involves active and reactive power balance as well as voltage magnitudes and angles [21]. Note that the two first models are linear programming models, i.e., these are convex models, while the third model is nonlinear and non-convex due to the products among trigonometric functions and voltages in the power balance equations [22].

In the literature, these models are employed to operate BESS and renewable sources in power systems from the point of view of the economic dispatch approach. Some of these works are presented as follows. The authors of [18] presented a unique nodal model for optimal dispatching batteries in AC microgrids. The resulting linear programming model is solved using the CPLEX solver, and the numerical results show the effect of the load and renewable generation behavior in the battery demeanor regarding state-of-charge and power injection/absorption. Berglund et al. [23] presented a linear-integer programming model to minimize the total cost of the facility, including the peak power demand and the energy costs, with the aim of extending battery lifespan. The proposed model uses a unique nodal representation and its solution is reached by using the GAMS optimization package. The authors of [2,20,24,25] presented linear programming and integer-linear programming models for installing and operating BESS in power systems using the DC equivalent of the AC power flow problem. The main aims of these works are the minimization of the investment costs in power systems regarding new transmission lines and the reduction of the greenhouse gas emissions produced by conventional thermal power sources. Most of these models have been solved using the GAMS software. The authors of [17,26,27] evaluated the complete AC multi-period optimal power flow model for the operation of BESS, considering a daily demand curve and renewable generation penetration including apparent power capabilities of the electronic convertors that are interconnected to the BESS. The authors of [21] proposed a nonlinear programming model with continuous an integer variables to install and operate batteries in distribution grids to reduce the total costs of the power losses during daily operation. The batteries are modeled with a binary behavior where the power injection is fixed regarding their maximums and minimums. The solution of this model is reached by implementing a Chu–Beasley genetic algorithm. The authors of [1] proposed an optimization approach to locate and operate BESS and distributed generators in medium- and low-voltage distribution networks, using a hybrid optimization approach that combines the simulating annealing method in the planning stage and a conic decomposition in the operation stage. Three test feeders composed of 11, 135, and 230 nodes are used to validate the proposed optimization methodology, considering the minimization of the investment, operation, and maintenance costs during the planning period. Numerical results demonstrate the effectiveness and robustness of the optimization approach regarding the planning period; however, the authors did not consider terms related to the environmental impact associated with greenhouse gas emissions. In addition, due to the usage of heuristics in the planning stage, the possibility of finding the global optimal is not ensured.

After the revision of the state-of-the-art associated with the optimal integration of BESSs systems in AC distribution grids, the main contributions of this study can be summarized as follows:

• ✓. The multi-objective formulation of the problem regarding the optimal operation of BESS in AC radial distribution networks using the branch optimal power flow representation, considering the simultaneous minimization of the CO${}_2$ gas emissions and the costs of the daily energy losses, is presented.

• ✓. The Pareto front is constructed using the multi-objective optimization approach via pondering factors by exploiting the potentialities of the GAMS software for nonlinear optimization.

• ✓. The different effects that voltage control in the substation and the active and reactive power injection in the BESSs have in the formation of the Pareto front were evaluated.

It is worth mentioning that the scope of this research is focused on the evaluation of different operative scenarios for multi-objective operation of BESSs in radial distribution networks, considering the linear representation of the state-of-charge in batteries and the amount of active power injected (or absorbed) to (or from) the grid, combined with the branch optimal power flow formulation, which has clearly not been proposed in the existing literature. However, it is important to note that even the optimal integration/operation of BESSs in electrical distribution grids can produce important technical and economical profits to the distribution system operators. Such an important fact regarding batteries (i.e., temperature-dependence and time-degradation) must be considered when batteries are installed in distribution grids, since the total expected benefit during its lifespan can be reduced significantly in comparison with ideal operative scenarios. This is an important area of research that the authors of [28] explored for direct current networks, and it can be extrapolated to AC grids in future works using our optimization model that serves as a starting point for future research.

The remainder of this document is organized as follows. Section 2 presents the complete multi-objective formulation of the problem of the optimal operation of BESSs in AC radial distribution networks by using the branch power flow formulation. Section 3 presents the solution methodology to solve the proposed multi-objective optimization problem by using the GAMS software that employs the pondering objective functions strategy. Section 4 shows the main characteristics of the 33-bus test feeder regarding its topology, demand behavior, renewable energy penetration, and BESSs’ locations, among others. Section 5 presents the computational implementation of the proposed multi-objective optimization model to operate BESSs in radial AC grids as well as the main results with their corresponding analysis and discussion. Section 6 presents the main concluding remarks derived from this research as well as possible future directions.

2. Multi-Objective Optimization Problem

The problem of the optimal operation of battery energy storage systems in AC grids is a nonlinear programming (NLP) problem. Here, we adopt the formulation of the power flow problem for strictly radial distribution networks that was proposed by Farivar and Low [19], by extending it to the multi-period domain and including the linear model of the batteries [29]. The complete formulation of the NLP problem is presented below.

2.1. Objective Functions

To operate a BESS in AC distribution systems, we consider two objective functions: the minimization of the CO${}_2$ greenhouse gas emissions and the daily cost of the energy losses in branches. These objective functions are defined as follows:

(1)$\begin{array}{c}\hfill min{z}_1=C{O}_2^{\mathrm{emissions}}\sum _{t\in \mathcal{T}}\sum _{i\in \mathcal{N}}{p}_{0i,t},\end{array}$

(2)$\begin{array}{c}\hfill min{z}_2=C_{ave}^{\mathrm{energy}}\sum _{t\in \mathcal{T}}\sum _{ij\in \mathcal{L}}{R}_{ij}{I}_{ij,t}^2,\end{array}$

where $z_1$ is the value of the objective function associated with the total emissions of CO${}_2$ to the atmosphere by diesel generators, $z_2$ is the cost of the total energy losses in all the branches of the energy, $C{O}_2^{\mathrm{emissions}}$ is the rate of CO${}_2$ emissions, $p_{0i,t}$ is the active power generation that leaves the substation, i.e., root node, to node i, $C_{ave}^{\mathrm{energy}}$ is the average cost of the energy, $R_{ij}$ is the resistance associated with the distribution line that connects nodes i and j, and $I_{ij,t}$ represents the magnitude of the current that flows through the line $ij$ at the period of time t. Observe that $\mathcal{N}$ and $\mathcal{T}$ are the sets that contain all the nodes of the grid and the periods of time of the study horizon.

2.2. Set of Constraints

The set of constraints associated with the operation of the BESS in electrical AC grids is composed of two main subsets, which are associated with the conventional power flow equations and the battery operation characteristics. Equations (3)–(6) present the conventional power flow equations formulated via branch model as presented in [30]:

(3)$\begin{array}{c}\hfill p_{ij,t}-R_{ij}{I}_{ij,t}^2-\sum _{k:\left(j,k\right)\in \mathcal{L}}{p}_{jk,t}=P_{j,t}-p_{j,t}^{dg}-p_{j,t}^b,\\ \hfill \left\{j\in \mathcal{N},t\in \mathcal{T}\right\}\end{array}$

(4)$\begin{array}{c}\hfill q_{ij}-X_{ij}{I}_{ij,t}^2-\sum _{k:\left(j,k\right)\in \mathcal{L}}{q}_{jk,t}=Q_{j,t}-q_{j,t}^{dg}-q_{j,t}^b,\\ \hfill \left\{j\in \mathcal{N},t\in \mathcal{T}\right\}\end{array}$

(5)$\begin{array}{c}\hfill V_{j,t}^2=V_{i,t}^2-2\left(R_{ij}{p}_{ij,t}+X_{ij}{q}_{ij,t}\right)+\left(R_{ij}^2+X_{ij}^2\right)I_{ij,t}^2,\\ \hfill \left\{(i,j)\in \mathcal{L},t\in \mathcal{T}\right\},\end{array}$

(6)$\begin{array}{c}\hfill I_{ij,t}^2=\frac{p_{ij,t}^2+q_{ij,t}^2}{V_{i,t}^2},\left\{(i,j)\in \mathcal{L},t\in \mathcal{T}\right\},\end{array}$

where $p_{ij,t}$ and $q_{ij,t}$ are the active and reactive power flows leaving node i to travel towards node j in a period of time t. $p_{jk,t}$ and $q_{jk,t}$ have the same definition applied to nodes j and k, respectively. Further, $P_{j,t}$ and $Q_{j,t}$ are the active and reactive power consumption measures at node j in a period of time t modeled as constant power loads. $p_{j,t}^{dg}$ and $q_{j,t}^{dg}$ are the active and reactive power generations in the distributed sources in the node j at the period of time t. Moreover, $p_{j,t}^b$ and $q_{j,t}^b$ correspond to the active and reactive power injections/absorptions in the battery b connected at node j in the period of time t and $V_{i,t}$ and $V_{j,t}$ are the voltage magnitudes in nodes i and j at the period of time t, respectively.

Remark 1.

Observe that the set of constraints in (3)–(6) is nonlinear and non-convex due to the presence of square variables associated with voltages, currents, and powers. However, the main complication of this model corresponds to Expression (6), since it defines the hyperbolic relation between voltages and powers [19].

The subset of constraints associated with batteries, element capacities, and voltage regulation bounds is described as follows:

(7)$\begin{array}{c}\hfill so{c}_{j,t}^b=so{c}_{j,t-1}^b-{\phi }_j^b{p}_{j,t}^b\Delta T,\left\{t\in \mathcal{T},j\in \mathcal{N}\right\},\end{array}$

(8)$\begin{array}{c}\hfill {\left(p_{j,t}^b\right)}^2+{\left(q_{j,t}^b\right)}^2\le {\left(s_j^{b,max}\right)}^2,\left\{t\in \mathcal{T},j\in \mathcal{N}\right\},\end{array}$

(9)$\begin{array}{c}\hfill p_{j,t}^{b,min}\le p_{j,t}^b\le p_{j,t}^{b,max},\left\{t\in \mathcal{T},j\in \mathcal{N}\right\},\end{array}$

(10)$\begin{array}{c}\hfill so{c}_j^{b,min}\le so{c}_{j,t}^b\le so{c}_j^{b,max},\left\{t\in \mathcal{T},j\in \mathcal{N}\right\},\end{array}$

(11)$\begin{array}{c}\hfill p_{j,t}^{gd,min}\le p_{j,t}^{gd}\le p_{j,t}^{gd,max},\left\{t\in \mathcal{T},j\in \mathcal{N}\right\},\end{array}$

(12)$\begin{array}{c}\hfill {\left(p_{j,t}^{gd}\right)}^2+{\left(q_{j,t}^{gd}\right)}^2\le {\left(s_j^{gd,max}\right)}^2,\left\{t\in \mathcal{T},j\in \mathcal{N}\right\},\end{array}$

(13)$\begin{array}{c}\hfill V_{j,t}^{min}\le V_{j,t}\le V_{j,t}^{max},\left\{t\in \mathcal{T},j\in \mathcal{N}\right\},\end{array}$

where $so{c}_{j,t}^b$ is the state-of-charge of the BESS connected at node j in the period of time t, and ${\phi }_j^b$ is the battery charging/discharging coefficient. Besides, $so{c}_j^{b,min}$ and $so{c}_j^{b,min}$ are the lower and upper limits of the state-of-charge variable (these are lower than one and greater than zero to increment battery lifespan; for Ion-Lithium batteries, these are recommended to be $10\%$ to $90\%$ [6]). Note that $p_{j,t}^{b,min}$ and $p_{j,t}^{b,max}$ are the active power limits of the batteries. As such, $s_j^{b,max}$ is the upper limit regarding power transference of the BESS installed at node j, $p_{j,t}^{gd,min}$ and $p_{j,t}^{gd,max}$ are the lower and upper power limits of the distributed generators (DGs) connected at the node j in the time period t, $s_j^{gd,max}$ is the maximum apparent power capability of the DG connected at node j, and $V_{j,t}^{min}$ and $V_{j,t}^{max}$ represent the minimum and maximum voltage regulation bounds admissible in all nodes of the network.

2.3. Interpretation of the Mathematical Model

The multi-objective optimization model defined in (1)–(13) has the following interpretation: Equations (1) and (2) represent the objective functions associated with the minimization of the CO${}_2$ greenhouse gas emissions, and the cost of the daily energy losses. Equations (3) and (4) correspond to the power balance equation at each node of the network, except in the voltage-controlled source. Equations (5) and (6) define the voltage drop at each branch and the hyperbolic relation between voltages and powers, respectively. Equation (7) determines the linear relationship between the state-of-charge and the power injection/absorption at each BESS [18]. Inequality constraints (8)–(10) define the apparent power circle of the battery, the maximum and minimum bounds associated with the active power injection/absorption, and state-of charge, respectively. Inequality expressions (11) and (12) define the distributed generation limits regarding apparent and active power generations. Finally, Expression (13) is known as the voltage regulation constraint, which defines the admissible voltage variations in all the nodes of the network.

Remark 2.

The solution of the mathematical model, (1)–(13), requires powerful optimization tools since this is a nonlinear, non-convex optimization problem with multiple variables defined in the continuous domain. Here, to solve this NLP problem, we adopt the implementation presented in [31], where the GAMS optimization package has been employed to solve a single-objective optimization problem regarding BESS selection and location. This implies that, in this paper, we extend the application of this software for multi-objective optimization problems as described in [2] for power system applications.

3. Solution Methodology

To resolve the nonlinear programming model defined in (1)–(13), the GAMS software was selected in this research. This software can solve different single-objective optimization problems starting from linear programming, nonlinear programming, and, in general, mixed-integer nonlinear programming problems [2]. However, we adapted this software via loops to solve the proposed NLP problem for multi-objective operation of batteries using the well-known pondering factors approach [32,33]. Table 1 summarizes some studies that worked with the GAMS software to solve different optimization problems to demonstrate its potential applications.

In Table 1, we can observe that, in the literature, there are multiple optimization problems solved with the GAMS software. This demonstrates that this is an excellent tool to address complex programming models such as the case of the multi-objective operation of BESSs in AC grids study in this research.

In general, the implementation of an optimization model in the software GAMS has the following main steps [50]:

• Define the sets associated with the groups of variables of the problem, i.e., set of periods of time $\mathcal{T}$, set of nodes $\mathcal{N}$, and set of branches $\mathcal{L}$.

• Define the scalars, parameters (vectors), and tables (matrices), i.e., active and reactive power demands ($P_{j,t}$ and $Q_{j,t}$), resistances and inductances per distribution line ($R_{ij}$ and $X_{ij}$), and the maximum and minimum bounds of the variables (i.e., $p_{j,t}^{b,min}$, $p_{j,t}^{b,max}$, $so{c}_j^{b,min}$, $so{c}_j^{b,max}$, and so on).

• Define the variables and their natures, i.e., continuous, binary, or discrete.

• Redact the equation names associated with each of the expressions in the optimization model, as well as the mathematical formulation of these equations using symbolic structure.

• Select the direction of the optimization, i.e., minimization, and the nature of the problem for being solved, i.e., NLP.

Remark 3.

The implementation of an optimization model in the software GAMS requires some familiarization with its syntax and general structure, since it works using plain text to formulate the complete mathematical model. For more details about the GAMS usage, see the works by Soroudi [2] and Castillo et al. [50].

A schematic implementation of an optimization model using the GAMS software is depicted using a general flow chart in Figure 1.

4. Information of the Test Feeders

This section presents the parametric information of the test feeders considered for validating the multi-objective operation of BESSs in rural distribution grids fed by diesel sources. The first test feeder corresponds to the IEEE 33-node test feeder and the second is the IEEE 69-node test feeder. All information regarding both distribution grids is presented below.

4.1. IEEE 33-Node Test Feeder

The IEEE 33-bus test feeder, which is a radial distribution network composed of 33 buses and 32 distribution lines, is operated at 12.66 kV at the substation bus (i.e., bus 1). The electrical connection between buses in the 33-bus test feeder is presented in Figure 2 and the parametric information of this test feeder is presented in Table 2. It is important to highlight that the information presented in this table corresponds to the load peak consumption in all the nodes of the system, which is surmised of 3715 kW and 2300 kvar. In addition, it is considered here that this medium-voltage distribution network is located in a rural area, where the main substation corresponds to a diesel generator, as reported in [27].

In this electrical network, the existence of four distributed generators is also considered: two of them arise from solar photovoltaic technology and two of them work with wind generation technology. The location of these distributed generators is as follows: (i) the PV${}_1$ source is installed at node 13 with a capacity of generation defined as 450 kW and the PV${}_2$ source is located at node 25 with a capability of 1500 kW; and (ii) the WT${}_1$ source is connected at the same node as the PV${}_1$ source with a nominal capacity of 825 kW and the PV${}_2$ is sited at node 30 with a rate of generation defined as 1200 kW. The information associated with the distributed sources and demand is listed in Table 3.

It is worth mentioning that the information regarding renewable generation presented in Table 3 was obtained by using an artificial neural network (ANN) to predict the daily generation considering a year of historic information that contains data about solar radiation, wind speed, humidity, pressure, and time. The complete information related to the implementation of the ANN can be found in [38], and some aspects regarding the ANN implementation are provided in Appendix A.

Regarding the location and sizing of the BESSs, we consider the information reported in [27] that includes three batteries allocated in this test feeder: (i) a BESS with an energy capacity of 1000 kWh that can be charged/discharged in 4 h (Battery Type A); (ii) a BESS with an energy capacity of 1500 kWh that can be charged/discharged in 4 h (Battery Type B); and (iii) a BESS with an energy capacity of 2000 kWh that can be charged/discharged in 5 h (Battery Type C). Observe that the power injection/absorption of the battery can be calculated by dividing the energy rate by the total charging/discharging time, i.e., for the Battery Type A, this power injection/absorption is 250 kW. In addition, to define the location of these batteries, we consider the solution reported in [31] that defines that Battery Types A–C are located at nodes 14, 31, and 6, respectively.

To evaluate the objective functions associated with the minimization of the greenhouse emissions and the daily energy loss costs, we consider that the diesel generator located at node 1 emits about $612.35$ kg/MWh and that the average energy costs are assumed to be US$ $0.1390$, as taken from [26], by converting Colombian pesos (COP$ $479.3389$) to American dollars.

4.2. IEEE 69-Node Test Feeder

The IEEE 69-bus test feeder is a radial distribution network composed of 69 buses and 69 distribution lines, which is operated at 12.66 kV at the substation bus (i.e., bus 1). The electrical connection between buses in the 69-bus test feeder is presented in Figure 3 and the parametric information of this test feeder is presented in Table 4. It is important to highlight that the information presented in this table corresponds to the load peak consumption in all nodes of the system, which surmises 3890.70 kW and 2693.6 kvar. In addition, here, it is considered that this medium-voltage distribution network is located in a rural area, where the main substation corresponds to a diesel generator and as reported in [27].

In this electrical network, the existence of four distributed generators is also considered, two of which work with solar photovoltaic technology and two of which work with wind generation technology. The location of these distributed generators is as follows: (i) the PV${}_1$ source is installed at node 12 with a capacity of generation defined as 1050 kW and the PV${}_2$ source is located at node 22 with a capability of 850 kW; and (ii) the WT${}_1$ is sited at node 61 with a rate of generation defined as 760 kW and the WT${}_2$ source is connected to the same node as the PV${}_1$ source with a nominal capacity of 1000 kW. Note that the information of the distributed sources and demand is listed in Table 3.

As this test feeder is considered an additional BESS regarding the batteries considered for the IEEE-33 node test feeder (i.e., Battery Type D), as reported in [27], this is a BESS with an energy capacity of 3000 kWh that can be charged/discharged in 6 h. The location of the BESS Types A–D is at nodes 40, 64, 16, and 9, respectively.

5. Numerical Simulations

This section shows the computational implementation of the proposed multi-objective optimization approach to operate BESSs in AC distribution systems considering the minimization of the greenhouse gas emissions versus costs of daily energy losses. We implemented this approach on a personal computer AMD Ryzen 7 3700U, 2.3 GHz, 16 GB RAM with 64-bits Windows 10 Home Single Language using the GAMS optimization environment.

The construction of the Pareto fronts presented in this section was made by using the weighting factors methodology that combines both objective functions as follows: $z_{\mathrm{comb}}^m=\vartheta z_1+(1-\vartheta )z_2$, where $\vartheta =m/20$ is the pondering factor. Note that m defines the number of combined objective function evaluations. In the proposed study, the parameter m begins from 0 in steps of 1–20. It is worth mentioning that the value of $\vartheta$ determines the number of solutions contained in the Pareto front.

For numerical purposes, different simulation cases in this section are presented and applied to the IEEE 33-bust test feeder and the IEEE 69-bus test feeder. For the first distribution grid, we considered the following cases: (i) effect of the substation voltage control; (ii) effect of the reactive power compensation with battery converters; and (iii) effect of the renewable energy variation in the Pareto front conformation. However, for the IEEE 69-node test feeder, to avoid repeating the analysis made in the IEEE 33-node test feeder, we only evaluated the operation of the batteries with unity and variable power factors, given the generators operate in their nominal condition.

5.1. IEEE-Bus Test Feeder

In this subsection, all the simulation scenarios considered for the IEEE-33 bus test feeder are presented, which include the effect of the substation voltage control, the variable power factor capabilities in batteries, and the possible variation of the renewable energy resources.

5.1.1. Effect of the Substation Voltage Control

To demonstrate the effect that has the assignation of the substation voltage in a predefined value, e.g., ${\mathbb{V}}_1=1.00\angle 0$, we simulated two cases. The first case considers that the substation is controlling the voltage magnitude at $1.00$ pu and the second case leaves this voltage magnitudes free, i.e., it can take values from $v_{min}$ to $v_{max}$. Figure 4 presents the Pareto fronts for both simulation cases.

The results in Figure 4 show the following. In the first case, the minimum and maximum values regarding the costs of energy losses are US$ $132.0450$ and US$ $157.1888$, which correspond to total CO${}_2$ emissions of about $13.8714$ and $6.5502$ tons/day, respectively. The difference between the extreme solutions are US$ $25.1438$ and $7.3212$ tons/day, respectively. In the second case, the minimum and maximum values regarding the costs of energy losses are US$ $108.1019$ and US$ $129.4134$, which correspond to total CO${}_2$ emissions of about $13.8294$ and $6.4853$ tons/day, respectively. The difference between the extreme solutions are US$ $20.6115$ and $7.3441$ tons/day, respectively. Finally, the main effect of maintaining the voltage magnitude of the substation in a fixed value (i.e., Case 1), in contrast to the solution reported by Case 2, is associated with an increment in the daily cost of the energy losses about US$ $26.4417$. However, when CO${}_2$ emissions are analyzed, we observe that the average difference is less than $0.0500$ tons/day. These results imply that fixing the voltage substation has a direct effect on the energy losses, and that this effect is negative, as shown in Figure 4.

The complete list of the objective function pairs is reported in Table 5, based on the information in Figure 4 from left to right.

Note that results in Table 5 demonstrate the multi-objective nature of the problem of the optimal operation of BESSs in AC grids, since the improvement of one of the objectives implies the deterioration of the other. In addition, observe that Solutions 6–12 can be considered attractive solutions to be selected by the distribution system operator, since these are in the intermediate area of the Pareto, which implies that both objectives perform well regarding their extremes. For example, Solution 9 in Case 1 produces a reduction of about $43.05\%$ regarding CO${}_2$ emissions and $7.11\%$ in relation to the maximum energy daily costs; in Case 2, these reductions are $43.08\%$ and $7.49\%$, respectively.

5.1.2. Effect of the Reactive Power Compensation with Battery Converters

In the previous subsection, we evaluate the effect of the voltage control in the substation bus considering unity power factor in the operation of the batteries and distributed generators. In this simulation case, we consider that the substation has free voltage magnitude and the batteries can support reactive power with their converters by assuming that the DGs will continue working with unity power factor [26]. In Figure 5, the Pareto front for the operational BESS, considering the variable power factor, is presented.

From the results in Figure 5, we can observe the following to be true. First, the extreme solutions of the Pareto front correspond to minimum CO${}_2$ emissions of $6.3136$ tons/day with a maximum daily energy costs of about US$ $59.1780$, while the maximum value of the CO${}_2$ emissions is $13.9245$ tons/day with minimum daily operative costs of US$ $37.5963$. Second, the difference between the extremes corresponds to a reduction of $54.66\%$ regarding CO${}_2$ emissions, while the costs of daily energy losses is about $34.47\%$. Third, when comparing the Pareto front in Figure 5 with Case 2 in the Pareto front presented in Figure 4, it is possible to observe that the injection of the reactive power in batteries allows the reduction of the daily operation costs by about US$ 70 per day, which represents a daily profit of $54.27\%$ to the grid operator.

To illustrate the behavior of the state-of-charge in batteries, we plot this variable for all the batteries in both extremes (see Figure 6) of the Pareto front in Figure 5.

The behaviors of the BESSs in Figure 6 show that, depending on the objective function prioritized, the BESSs will work with different charging/discharging profiles along the day. Depending on the minimization interest, they decide to store or provide energy in different periods of time (see comparison of batteries in nodes 6 and 14), especially in the period of time between 14 and 28. In addition, note that, in the case of minimizing the daily energy costs, the BESSs start to provide energy continuously from the period of time 16 to the period of time 46, while, in the case of the CO${}_2$ minimization, these behaviors do not occur. This is attributable to the availability of the renewable resources and the demand behavior, which has two main peaks that require additional power injections from the BESS to minimize the amount of power injection in the conventional source that is directly connected with the amount of CO${}_2$ emissions, as seen in Equation (1).

5.1.3. Effect of the Renewable Energy Variation in the Pareto Front Conformation

To present the effect of the renewable energy variation due to the stochastic nature of the primary energy resources, i.e., wind speed and solar radiation, we now present a simulation case where the amount of active power injection from the renewable energy resources varies from $50\%$ to $100\%$ for each varying level of renewable energy availability. As such, the Pareto front presented in Figure 7 is obtained, considering unity power factor in the renewable generation, active, and reactive power injections in the BESS, and the free behavior in the slack voltage magnitude.

The following can be noted from the behavior of the Pareto fronts in Figure 6. First, depending on the availability of the renewable energy resources during a day of operation, the Pareto front can be very restricted. For instance, in the case of a renewable energy availability of $50\%$, the daily operative costs are between US$ 55 and US$ 58, while CO${}_2$ emissions are between 20 and 21 tons/day. Second, when the amount of active power injection with renewable energy resources increases, the behavior of the Pareto front regarding the CO${}_2$ emissions presents higher differences in their extremes. For example, for an availability of $50\%$, the difference between the extreme solutions is about $1.3506$ tons/day, and, for a $90\%$ renewable generation availability, this value is about $2.9869$ tons/day. This implies the growing importance of decisions-making regarding the operation of the distribution network on the part of the distribution company when the amount power injection from renewable generation increases. Finally, the biggest effect associated with the availability of the renewable generation to inject active power corresponds to the amount CO${}_2$ emissions, since this variable is directly connected to the total power injection in the slack node, which is then directly connected to the power injections in the distributed generators. This is particularly true while the daily operative costs for all the percentages of penetration between the extremes remain the same in the case of the nominal operation, i.e., 100% of renewable generation availability. These results imply that the BESS can control the grid losses in better form, since it depends on the voltage drops in lines, as compared to the the case of the greenhouse emissions, since these depend directly on the usage of the diesel generator.

5.2. IEEE 69-Bus Test Feeder

For the IEEE 69-bus test feeder, we consider the possibility of operating the BESSs while considering that the voltage control of the substation is free and can can be moved from 0.90 to 1.10 pu by maintaining its angle at zero for all periods of time. The Pareto front for this test feeder is presented in Figure 8. Note that the extreme solutions in the Pareto front are the following: (i) the minimum daily energy costs is US$ 65.8980 with maximum CO${}_2$ emissions of 18.5373 tons/day; and (ii) the minimum emissions of CO${}_2$ are 10.9535 tons/day, with a maximum daily costs of energy losses of US$ 89.1888.

It is worth mentioning that the difference between the extreme solutions in the Pareto front in Figure 7 is US$ 23.2908 and 7.5838 tons/day. These values imply that, depending on the distribution system operator, the daily indicators can reduce until these values are selected as a function of the operative point. In addition, it is possible to observe that the most promissory solutions that allow minimizing both objective functions are in the region between 12 and 15 tons/day, since these have daily energy loss costs of between US$ 79 and US$ 71.

6. Conclusions

We studied the problem of operation of a BESS in radial AC grids from the multi-objective point of view, considering the simultaneous minimization of CO${}_2$ emissions and the costs of the daily energy losses. The GAMS software combined with the pondering objective functions approach were employed to solve the resulting nonlinear, non-convex, and multi-period optimal power flow problem. Numerical validations demonstrated that the assignation of the voltage magnitude in the substation bus deteriorates the behavior of the optimal Pareto front since both objectives present higher values in comparison with the scenario where voltage magnitude in the substation is left free.

When reactive power injections are available for operating the BESS, considering its power electronic converters, the objective function regarding energy loss costs is drastically reduced. This was made possible since these reactive power injections allowed the improvement of voltage profiles, which is directly connected with the reduction of magnitude of the current in branches. This reduction was about $54.27\%$ in comparison to the unity factor power case.

The effect of the renewable generation availability was tested for different percentages of active power generation in distributed sources. This demonstrated that the amount of greenhouse gas emissions is directly connected to the power supply by distributed generators, while daily energy losses can be effectively controlled for the BESSs by maintaining these on the extreme values of the nominal operative case, i.e., between US$ $37.5963$ and US$ $59.1780$, for all the renewable generation penetration cases evaluated.

In the IEEE 69-bus system, it was observed that the maximum reduction of CO${}_2$ emissions is about 40.91% when the daily cost of energy losses has a maximum value, i.e., US$ 89.1888. In the case of energy losses, the maximum reduction is about 26.11% when the maximum CO${}_2$ emissions occur. These values correspond to the extreme solutions in the Pareto front in Figure 7, which implies that, depending on the distribution system operator’s decision, the maximum improvements possible will be contained in the aforementioned percentages of reduction.

The following research can be developed as part of future research efforts following this paper’s contribution. First, the reformulation of the proposed branch multi-period optimal power flow model for multi-objective operation of the BESS in AC grids with radial structure into a convex conic representation that ensures the global optimum finding from the formal mathematical point of view. Second, the application of proposed multi-objective optimization model to direct current networks with high penetration of renewable energy resources and batteries. Finally, the inclusion of new battery behaviors in the proposed optimization model, such as lifespan and deep discharging characteristics, among others, to approximate the best behavior of these in real operative conditions.

Author Contributions

Conceptualization, F.M.-M., O.D.M., L.F.G.-N., J.C.H., and C.A.R.-V.; Methodology, F.M.-M., O.D.M., L.F.G.-N., J.C.H., and C.A.R.-V.; Investigation, F.M.-M., O.D.M., L.F.G.-N., J.C.H., and C.A.R.-V.; and Review—writing and editing, F.M.-M., O.D.M., L.F.G.-N., J.C.H., and C.A.R.-V. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Agencia Estatal de Investigación, Spain (AEI) and the Fondo Europeo de Desarrollo Regional (FEDER) aimed at the Challenges of Society (Grant No. ENE 2017-83860-R “Nuevosservicios de red para microredes renovables inteligentes. Contribución a la generación distribuida residencial”).

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments

This work was supported in part by the Centro de Investigación y Desarrollo Científico de la Universidad Distrital Francisco José de Caldas under Grant 1643-12-2020 associated with the project, “Desarrollo de una metodología de optimización para la gestión óptima de recursos energéticos distribuidos en redes de distribución de energía eléctrica”, and in part by the Dirección de Investigaciones de la Universidad Tecnológica de Bolívar under Grant PS2020002 associated with the project, “Ubicación óptima de bancos de capacitores de paso fijo en redes eléctricas de distribución para reducción de costos y pérdidas de energía: Aplicación de métodos exactos y metaheurísticos”.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Renewable Generation Forecasting

The renewable generation dependence on the weather conditions associated with solar radiation and temperature for photovoltaic plants make this source non-dispatchable, which clearly affects the grid operation, causing increments in the operative costs, low voltage profile, and high energy losses [26]. For this reason, to operate electric distribution networks in rural and urban areas, it is necessary to have an adequate short-term forecasting approach to estimate the renewable energy output for the next day. In the literature, some approaches are proposed to predict the renewable energy behavior. These include fuzzy logic [52], stochastic optimization methods [53], similar day approach [54], and artificial neural networks [26,55], among others. Here, we employed the methodology proposed in [26] to forecast the renewable energy for the day-ahead operation of the distribution grid.

It is worth mentioning that, for the sake of simplicity, in the ANN implementation, we only present information of the PV prediction. However, the complete details for wind generation forecasting can be seen in [55].

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. General GAMS implementation for solving NLP optimization problems.

Figure 2. IEEE 33-bus test feeder configuration.

Figure 3. IEEE 69-bus test feeder configuration.

Figure 4. Behavior of the greenhouse gas emissions vs. daily energy loss costs depending on the voltage magnitude at the substation node.

Figure 5. Behavior of the greenhouse gas emissions vs. daily energy loss costs considering reactive power injections in the BESS.

Figure 6. Behavior of the state-of-charge in the BESS for the extreme solutions of the Pareto front in Figure 5: (a) minimum CO2 emissions; and (b) minimum daily energy costs.

Figure 7. Performance of the Pareto front for different levels of renewable generation availability.

Figure 8. Performance of the Pareto front for the IEEE 69-bus system.

Table 1

Engineering optimization problems solved in the software GAMS.

Optimization Problem	References
Optimal location and sizing distributed generation in AC grids	[34,35,36,37,38]
Distribution system planning	[39,40,41]
Optimal location of capacitor banks in distribution networks	[42,43,44]
Optimal location and operation of battery energy storage systems in distribution networks	[2,17,26,31]
Efficient design of osmotic generation plants	[45,46,47]
Economic dispatch of thermal plants in power systems	[2,48]
Solution of the general engineering problems using GAMS	[49,50,51]

Table 2

IEEE 33-bus test feeder information regarding branches and loads.

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

Table 3

Behavior of the renewable generation and demand for a daily operation.

Time (s)	PV${}_1$ (p.u)	PV${}_2$ (p.u)	WT${}_1$ (p.u)	WT${}_2$ (p.u)	Demand (p.u)
0.0	0	0	0.633118295	0.489955551	0.34
0.5	0	0	0.629764678	0.467954207	0.28
1.0	0	0	0.607259323	0.449443905	0.22
1.5	0	0	0.609254545	0.435019277	0.22
2.0	0	0	0.605557422	0.437220792	0.22
2.5	0	0	0.630055346	0.437621534	0.20
3.0	0	0	0.684246423	0.450949300	0.18
3.5	0	0	0.758357805	0.453259348	0.18
4.0	0	0	0.783719339	0.469610539	0.18
4.5	0	0	0.815243582	0.480546213	0.20
5.0	0	0	0.790557706	0.501783479	0.22
5.5	0	0	0.738679217	0.527600299	0.26
6.0	0	0	0.744958950	0.586555316	0.28
6.5	0	0	0.718989730	0.652552760	0.34
7.0	0.039123365	0.026135642	0.769603567	0.697699990	0.40
7.5	0.045414292	0.051715061	0.822376817	0.774442755	0.50
8.0	0.065587179	0.110148398	0.826492212	0.820205405	0.62
8.5	0.132615282	0.263094042	0.848620129	0.871057775	0.68
9.0	0.236870796	0.431175761	0.876523598	0.876973635	0.72
9.5	0.410356256	0.594273035	0.904128455	0.877065236	0.78
10.0	0.455017818	0.730402039	0.931213527	0.897955131	0.84
10.5	0.542364455	0.830347309	0.955557477	0.903245007	0.86
11.0	0.726440265	0.875407050	0.965504834	0.916903429	0.90
11.5	0.885104984	0.898815348	0.971037333	0.924757605	0.92
12.0	0.924486326	0.975683083	0.972218577	0.942224932	0.94
12.5	1	1	0.980049847	0.949956724	0.94
13.0	0.982041153	0.978264398	0.981135531	0.963773634	0.90
13.5	0.913674689	0.790055240	0.988644844	0.974977461	0.84
14.0	0.829407079	0.882557147	0.991393173	0.986750539	0.86
14.5	0.691912077	0.603658738	0.998815517	0.995058133	0.90
15.0	0.733063295	0.606324907	1	1	0.90
15.5	0.598435064	0.357393267	0.996070963	0.998107341	0.90
16.0	0.501133849	0.328035635	0.987258076	0.997690423	0.90
16.5	0.299821403	0.142423488	0.976519817	0.993076899	0.90
17.0	0.177117518	0.142023463	0.929542167	0.982629597	0.90
17.5	0.062736095	0.072956701	0.876413965	0.972084487	0.90
18.0	0	0.019081590	0.791155379	0.930225756	0.86
18.5	0	0.008339287	0.691292162	0.891253999	0.84
19.0	0.000333920	0	0.708839248	0.781950905	0.92
19.5	0	0	0.724074349	0.660094138	1.00
20.0	0	0	0.712881960	0.682715246	0.98
20.5	0	0	0.733954043	0.686617947	0.94
21.0	0	0	0.719897641	0.681865563	0.90
21.5	0	0	0.705502389	0.717315757	0.84
22.0	0	0	0.703007456	0.718080346	0.76
22.5	0	0	0.686551618	0.726890145	0.68
23.0	0	0	0.687238555	0.734452193	0.58
23.5	0	0	0.682569771	0.739699146	0.50

Table 4

IEEE 69-bus test feeder information regarding branches and loads.

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.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.3811	0.1941	40.4	30	40	41	0.7283	0.8509	102	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.0478	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	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	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.1140	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.3463	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	10	10	66	67	0.0047	0.0014	18	13
33	34	1.7080	0.5646	14	14	12	68	0.7394	0.2444	28	20
34	35	1.4740	0.4873	4	4	68	69	0.0047	0.0016	28	20

Table 5

Objective function values of the Pareto front in Figure 4.

Sol. No.	Case 1		Case 2
	CO${}_{\mathbf{2}}$ (Tons/day)	Losses (US$)	CO${}_{\mathbf{2}}$ (Tons/day)	Losses (US$)
1	13.8713	132.0450	13.8293	108.1019
2	11.6074	133.7900	11.5511	109.5534
3	10.6953	135.5350	10.632	111.0050
4	10.0223	137.2800	9.9557	112.4566
5	9.4704	139.0250	9.4108	113.9082
6	8.9971	140.7700	8.9491	115.3597
7	8.5828	142.5150	8.5438	116.8113
8	8.2202	144.2600	8.1877	118.2629
9	7.8995	146.0050	7.8712	119.7145
10	7.6145	147.7500	7.5868	121.1660
11	7.3612	149.4950	7.3313	122.6176
12	7.1359	151.2400	7.1025	124.0692
13	6.9356	152.9850	6.9008	125.5208
14	6.7571	154.7299	6.7216	126.9723
15	6.5982	156.4749	6.5651	128.4239
16	6.5501	157.1888	6.4852	129.4134