Optimal Economic–Environmental Operation of BESS in AC Distribution Systems: A Convex Multi-Objective Formulation

Abstract

This paper deals with the multi-objective operation of battery energy storage systems (BESS) in AC distribution systems using a convex reformulation. The objective functions are CO₂ emissions, and the costs of the daily energy losses are considered. The conventional non-linear nonconvex branch multi-period optimal power flow model is reformulated with a second-order cone programming (SOCP) model, which ensures finding the global optimum for each point present in the Pareto front. The weighting factors methodology is used to convert the multi-objective model into a convex single-objective model, which allows for finding the optimal Pareto front using an iterative search. Two operational scenarios regarding BESS are considered: (i) a unity power factor operation and (ii) a variable power factor operation. The numerical results demonstrate that including the reactive power capabilities in BESS reduces 200 kg of CO₂ emissions and USD 80 per day of operation. All of the numerical validations were developed in MATLAB 2020b with the CVX tool and the SEDUMI and SDPT3 solvers.

1. Introduction

Energy storage systems play a key role in technological and industrial development, with applications in the smartphone, computers, conventional and electric vehicles, distribution networks, and renewable electricity generation [1,2,3], among others. In general, these energy storage systems are based on chemical technologies; most commonly lead–acid [4], nickel–cadmium [5], and lithium-ion [6,7] batteries. In electrical grids, battery energy storage systems (BESS) are widely used for power and energy compensation [8], power oscillation damping [9], frequency generation, and voltage support [10]. One of the fundamental applications of BESS is the use of energy from renewable sources to reduce the operating costs of electrical grids and to mitigate greenhouse gas emissions [11,12].

Including BESS in electricity distribution networks is a complex task in terms of control and optimization. For control, efficient energy management strategies must be designed considering non-linear characteristics of the converters that interface them with the electrical grid, which requires applying non-linear control methods for an adequate operation of BESS and their local and or global dynamic stability [13,14]. Some of the control methods that are applied to BESS systems include sliding mode control [15], passivity-based control [16], exact feedback linearization [17], fuzzy logic [18,19], neural networks [20], and classic proportional–integral–derivative (PID) controllers [21], among others. For optimization, a major challenge to overcome with BESS is the economic and environmental dispatch, which is a non-linear and nonconvex optimization problem with a difficult solution. Generally, this problem requires determining BESS charging/discharging periods taking advantage of available renewable sources during a defined period of operation (typically 24 h) in order to reduce one or more objective functions. The main complicating factor of this problem is the presence of power balance equations in each node of the system, because these equations are non-linear and typically require numerical solution methods.

In line with the above, the authors of [22] proposed a non-linear programming model for operating batteries on alternating current (AC) electrical grids with a power factor of 1.0. This model was solved in a general algebraic modeling system (GAMS). An improvement of this model was proposed in [8], which considers the possibility of injecting reactive power; however, it was solved using GAMS software, which, given the non-linear nature of the problem, does not necessarily lead to the global optimum. The authors [23] proposed a solution technique for the battery location and operative problem in distribution networks using a genetic algorithm combined with a particle swarm algorithm. In this model, BESS was run with a unity power factor, and their power injection was modeled with binary variables, complicating their solution. The authors [11] presented a mixed-integer non-linear programming model to determine the battery location and size in AC networks, considering unity and variable power factors. The solution to this mathematical model was found in the GAMS software. Given the non-linear and nonconvex nature of the problem, the global optimum is not necessarily found. The authors [24] proposed a reformulation of the second-order cone programming model to operate batteries in AC networks, which guarantees finding the global optimum; however, the reformulation of the problem is based on a single-objective equivalent, combining greenhouse gas emissions with the cost of energy losses, which prevents a Pareto optimal front. The authors of [25] proposed a semidefinite programming (SDP) model for solving the problem of volt/VAR optimization in unbalanced distribution grids with BESS to minimize power loss. However, the problem of the SDP model lies in recovering the rank-one solution from the matrix, which complicates it at the time of a real application. The authors of [26] presented a methodology based on the equilibrium optimization algorithm for the optimal integration of photovoltaic (PV) with BESS in order to improve the distribution network performance. The authors of [27] showed an improved sunflower optimization algorithm for optimal capacitor banks re-allocation in smart distribution networks. The authors of [28] performed a methodology using the coyote optimization algorithm to reduce the losses in the electricity distribution network by employing the BESS. Other studies have proposed optimization models for series-parallel configurations of battery cells to balance their thermodynamic performance [29]; however, these models overlook the distribution network in their formulation and primarily focus on electric vehicle applications [21]. Additionally, it is also important to determine the parameters adequately for the BESS in order to apply some methods mentioned above; the authors of [30] described an enhanced procedure to estimate the BESS parameters using an equilibrium optimization algorithm.

Unlike previous studies, this article contributes the following:

✓

A multi-purpose convex formulation for the efficient operation of energy storage systems in alternating current distribution networks, considering the simultaneous minimization of carbon dioxide emissions in diesel plants and energy losses in conductors. To this end, a reformulation of the second-order cone programming model of the multi-period optimal flow problem is proposed, with the main advantage of ensuring that the global optimum is found;

✓

The application of the weighted weights method to transform the multi-objective optimization model into a single-objective equivalent, which is iteratively solved to obtain the Pareto optimal front by solving this model using the Self-Dual-Minimization (SeDuMi) and SemiDefinite Programming (SDPT3) tools available in the MATLAB CVX tool;

✓

The assessment of the effect of manageable reactive power injection in the batteries concerning the cost of CO₂ losses and emissions when compared with the unity power factor scenario.

The remainder of this document is organized as follows: In Section 2, the formulation of the multi-objective operation problem of BESS in distribution systems is presented, highlighting their non-linear and nonconvex characteristics given by power balance equations and by voltage drops in branches. In Section 3, a convex reformulation is proposed by transforming the original non-linear model into a second-order cone programming model. Section 4 briefly describes the solution method; this method combines objective functions using weighted weights. In Section 5, the main characteristics of the radial IEEE 33-node test system are presented, composed of a slack node, four distributed generators, and three batteries. In Section 6, the computational validation is presented, considering two operational scenarios for BESS, based on unity and variable power factors. In Section 7, the main conclusions of this study are presented, as well as future perspectives.

2. Mathematical Formulation

The problem of the optimal operation of battery energy storage systems in distribution networks is formulated as a multi-objective optimization problem in which mitigating the emissions of greenhouse gases (CO₂ emissions) and reducing the cost of energy losses during the operating period are the objective functions. The structure of these objective functions is presented below:

$$\begin{array}{c}\hfill min{f}_1=C{O}_2^{\mathrm{em}}\sum _{t\in \mathcal{T}}\sum _{i\in \mathcal{N}}{p}_{0i,t}\Delta t,\end{array}$$

(1)

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

(2)

where $f_1$ represents the value of the objective function associated with CO₂ emissions to the atmosphere produced by diesel generators (measured in CO₂ tons per day, i.e., tons/day) and $f_2$ is the value of the objective function related to the daily cost of energy losses in all conductors of the distribution network (measured in USD/day); ${\mathit{CO}}_2^{\mathrm{em}}$ represents the rate of carbon dioxide emissions for medium-capacity diesel generation (i.e., <10 MW); $\Delta t$ is the fraction of time associated with the study period, i.e., 30 min for this study; $C_{\mathrm{ave}}^{\mathrm{energy}}$ corresponds to the average kilowatt-hour cost of grid losses during electric power transmission; $p_{0i,t}$ represents the amount of active power flowing from the substation node (i.e., source node 0) to the node i in the time period t; $R_{ij}$ corresponds to the value of the resistance between the nodes i and j; $I_{ij,t}$ corresponds to the magnitude of the current that flows between the nodes i and j in the time period t. It should be noted that $\mathcal{N}$, $\mathcal{L}$, and $\mathcal{T}$ correspond to the sets that contain all system nodes, all network sections, and all time periods, respectively.

The set of constraints associated with the optimal multi-objective dispatch problem of BESS contains the classic equations of active and reactive power balance and voltage drops in branches [31], as well as element capacity constraints, in addition to battery operating characteristics. The complete set of constraints of the optimization problem under study is presented below.

$$\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 \left\{j\in \mathcal{N}, t\in \mathcal{T}\right\},\end{array}$$

(3)

$$\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 \left\{j\in \mathcal{N}, t\in \mathcal{T}\right\},\end{array}$$

(4)

$$\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 \left\{(i,j)\in \mathcal{L}, t\in \mathcal{T}\right\},\end{array}$$

(5)

$$\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}$$

(6)

$$\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}$$

(7)

$$\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}$$

(8)

$$\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}$$

(9)

$$\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}$$

(10)

$$\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}$$

(11)

$$\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}$$

(12)

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

(13)

where $p_{ij,t}$ and $q_{ij,t}$ represent the flows from the node i to the node j in the time period t; $p_{jk,t}$ and $q_{jk,t}$ are defined as the nodes j and k above, respectively; $P_{j,t}$ and $Q_{j,t}$ correspond to the active and reactive power consumptions at node j in each time period t. These loads have been modeled as constant power loads: $p_{j,t}^{dg}$ and $q_{j,t}^{dg}$ are the active and reactive powers generated in the distributed sources connected to the node j in the time period t; $p_{j,t}^b$ and $q_{j,t}^b$ represent the active and reactive power generation or consumption in the battery b connected to the node j in the time period t; $V_{i,t}$ and $V_{j,t}$ correspond to the voltage magnitudes at the nodes i and j in the time period t, respectively. $so{c}_{j,t}^b$ represents the state of charge of the battery connected to the node j in the time period t, with ${\phi }_j^b$ being the battery charge/discharge coefficient. $so{c}_j^{b,min}$ and $so{c}_j^{b,min}$ correspond to the lower and upper bounds for the state-of-charge battery variables (these values should be higher than zero and lower than one; more specifically, for lithium-ion batteries, the values ranging from 10% to 90% are recommended, according to [8]). $p_{j,t}^{b,min}$ and $p_{j,t}^{b,max}$ represent the active power bounds of the batteries, $s_j^{b,max}$ being the maximum power transfer/absorption bounds of the battery b connected to the node j. $p_{j,t}^{gd,min}$ and $p_{j,t}^{gd,max}$ represent the active and reactive power bounds of the renewable generation sources connected to the node j in the time period t, $s_j^{gd,max}$ being the maximum apparent power transfer, and $V^{min}$ and $V^{max}$ correspond to minima and maxima associated with voltage regulation in all nodes of the distribution network.

In the mathematical model defined from (1) to (13), Equations (3) to (6) are nonlinear and non-convex due to the presence of squared current, voltage, and power variables. However, as stated by [31], the main complication of this mathematical model is the constraint (6), which is associated with the hyperbolic relationship that exists between voltages, currents, and powers.

The mathematical model (1)–(13) that represents the multi-objective operation of energy storage systems in alternating current distribution networks is interpreted as follows: Equations (1) and (2) define the objective functions of interest, the first of which being related to the mitigation of greenhouse gas emissions associated with carbon dioxide, and the second being related to the daily cost of energy losses in conductors of the distribution network. Equations (3) and (4) define the balance between active and reactive power in each node of the system, whereas Equation (5) is related to the voltage drop in each branch of the distribution system. Equation (6) defines the non-linear relationship between voltage and current, according to Tellegen’s theorem, associated with the definition of electric power; Equations (7) to (10) define, in their order, the linear relationship between the state of charge and the active power injected/imported or absorbed with the battery, the apparent power transfer capacity, and the active power injection bounds and the admissible state-of-charge bounds, respectively. Equations (11) and (12) are associated with the active power generation capacity of renewable generators, as well as their apparent power transfer bounds, respectively. Lastly, inequality (13) is known as the voltage regulation constraint, which is responsible for limiting the maximum voltage drop of the network, according to regulatory requirements applicable to medium voltage levels.

A key characteristic of the mathematical optimization model proposed for the problem of the efficient operation of energy storage systems in energy distribution networks is the ability to use the converter that integrates the battery as a reactive power compensation element [32,33]. It is worth mentioning that reactive power compensation with the battery compensator is entirely plausible because this device essentially operates as a static power compensator, whereby, based on an efficient control method, reactive power can be managed in two quadrants as a function of electrical grid requirements [34]. However, the stage of the voltage source converter control for reactive power management is beyond the scope of this article, so readers should refer to the references [8,35] for further details.

Deriving an equivalent convex mathematical model of the multi-objective problem defined from (1) to (13) only requires finding the convex equivalent of Equations (3) to (6) because the other equations, including the objective functions, are convex. The quadratic objective function associated with the power losses is convex because the resistance parameters are always higher than zero [24].

The convexification of the mathematical model (1)–(13) will be described in detail in the following section using second-order cone programming, as discussed in [31].

3. Convex Reformulation

To ensure that the global optimum is found in the nonlinear formulation of (1) presented by (13) for the economic–environmental dispatch problem of batteries in AC distribution systems with radial configuration, in this study, this model is transformed into an equivalent second-order cone programming model. This reformulation of the problem basically consists of transforming the nonlinear Equations (3) to (6) into equivalent planes, i.e., affine constraints. For this purpose, the following auxiliary variables are considered: $l_{ij,i}=I_{ij,t}^2$ and $u_{j,t}=V_{j,t}^2$. With these definitions, the objective function associated with energy losses (see Equation (2)), active and reactive power balance constraints, i.e., Equations (3) and (4), become linear, as does voltage drop constraints, i.e., (5). Therefore, the only complicated constraint is Equation (6) because it contains a product between variables, as shown below:

$$\begin{array}{c}\hfill p_{ij,t}^2+q_{ij,t}^2=l_{ij,t}{u}_{j,t}, \left\{(i,j)\in \mathcal{L}, t\in \mathcal{T}\right\}.\end{array}$$

(14)

This restriction can be easily transformed using the equivalent hyperbolic relation of the product of variables, as detailed below (set notation is omitted for simplicity):

$$\begin{array}{cc}\hfill l_{ij,t}{u}_{j,t}=& \frac{1}{4}{\left(l_{ij,t}+u_{j,t}\right)}^2-\frac{1}{4}{\left(l_{ij,t}-u_{j,t}\right)}^2,\hfill \\ \hfill p_{ij,t}^2+q_{ij,t}^2=& \frac{1}{4}{\left(l_{ij,t}+u_{j,t}\right)}^2-\frac{1}{4}{\left(l_{ij,t}-u_{j,t}\right)}^2\hfill \\ \hfill {\left(l_{ij,t}+u_{j,t}\right)}^2=& {\left(2p_{ij,t}\right)}^2+{\left(2q_{ij,t}\right)}^2+{\left(l_{ij,t}-u_{j,t}\right)}^2,\hfill \\ \hfill ∥\begin{array}{c}2p_{ij,t}\\ 2q_{ij,t}\\ l_{ij,t}-u_{j,t}\end{array}∥=& l_{ij,t}+u_{j,t}, \left\{(i,j)\in \mathcal{L}, t\in \mathcal{T}\right\}.\hfill \end{array}$$

(15)

The constraint (15) remains nonconvex [24] due to the condition of equality; accordingly, to obtain a conic equivalent, this constraint is relaxed using the inequality symbol, as follows:

$$\begin{array}{c}\hfill ∥\begin{array}{c}2p_{ij,t}\\ 2q_{ij,t}\\ l_{ij,t}-u_{j,t}\end{array}∥\le l_{ij,t}+u_{j,t}, \left\{(i,j)\in \mathcal{L}, t\in \mathcal{T}\right\}.\end{array}$$

(16)

Note that, with Equation (16), the multi-objective mathematical model (1)–(13) becomes convex. This optimization is known as second-order cone programming (SOCP). The complete mathematical model is presented below, describing the constraints (8) and (12), as in (16), to standardize the structure of the mathematical model.

$$\begin{array}{c}\hfill \begin{array}{c}\hfill min{f}_1=C{O}_2^{\mathrm{em}}\sum _{t\in \mathcal{T}}\sum _{i\in \mathcal{N}}{p}_{0i,t},\\ \hfill min{f}_2=C_{\mathrm{ave}}^{\mathrm{energy}}\sum _{t\in \mathcal{T}}\sum _{ij\in \mathcal{L}}{R}_{ij}{l}_{ij,t}\Delta t,\\ \hfill p_{ij,t}-R_{ij}{l}_{ij,t}-\sum _{k:\left(j,k\right)\in \mathcal{L}}{p}_{jk,t}=P_{j,t}-p_{j,t}^{dg}-p_{j,t}^b, \left\{j\in \mathcal{N}, t\in \mathcal{T}\right\},\\ \hfill q_{ij}-X_{ij}{l}_{ij,t}-\sum _{k:\left(j,k\right)\in \mathcal{L}}{q}_{jk,t}=Q_{j,t}-q_{j,t}^{dg}-q_{j,t}^b, \left\{j\in \mathcal{N}, t\in \mathcal{T}\right\},\\ \hfill u_{j,t}=u_{i,t}-2\left(R_{ij}{p}_{ij,t}+X_{ij}{q}_{ij,t}\right)+\left(R_{ij}^2+X_{ij}^2\right)l_{ij,t}, \left\{(i,j)\in \mathcal{L}, t\in \mathcal{T}\right\},\\ \hfill ∥\begin{array}{c}2p_{ij,t}\\ 2q_{ij,t}\\ l_{ij,t}-u_{j,t}\end{array}∥\le l_{ij,t}+u_{j,t}, \left\{(i,j)\in \mathcal{L}, t\in \mathcal{T}\right\},\\ \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\},\\ \hfill ∥\begin{array}{c}{p}_{j,t}^b\\ q_{j,t}^b\end{array}∥\le s_j^{b,max}, \left\{t\in \mathcal{T}, j\in \mathcal{N}\right\},\\ \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\},\\ \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\},\\ \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\},\\ \hfill ∥\begin{array}{c}{p}_{j,t}^{gd}\\ q_{j,t}^{gd}\end{array}∥\le s_j^{gd,max}, \left\{t\in \mathcal{T}, j\in \mathcal{N}\right\},\\ \hfill {\left(V^{min}\right)}^2\le u_{j,t}\le {\left(V^{max}\right)}^2, \left\{t\in \mathcal{T}, j\in \mathcal{N}\right\},\end{array}\end{array}$$

(17)

The mathematical model (17) can be solved using the interior-point or logarithmic barrier method, which, according to [31,36], guarantees that the global optimum is found thanks to the convex structure of the entire mathematical model.

To ensure that the initial and final states of the batteries reach reference values, the following battery life constraints are added to the mathematical model (17):

$$\begin{array}{c}\hfill so{c}_{j,t=t_i}^b=so{c}_j^{b,\mathrm{initial}}, \left\{t_i\in \mathcal{T}, j\in \mathcal{N}\right\},\\ \hfill so{c}_{j,t=t_f}^b=so{c}_j^{b,\mathrm{final}}, \left\{t_f\in \mathcal{T}, j\in \mathcal{N}\right\},\end{array}$$

where $t_i$ and $t_f$ are the start and end times of the study period, and $so{c}_j^{b,\mathrm{initial}}$ and $so{c}_j^{b,\mathrm{final}}$ represent the initial and final state of charge of the battery b connected to the node j. In this study, these values are assumed to be 50%, as recommended in [8].

To classify and characterize the convex optimization model (17), we present a complete classification table where the number of variables and their nature and the type and number of constraints, including the objective function and its components, are listed. Note that n (number of nodes) is defined as the cardinality of the set $\mathcal{N}$, and p (number of periods) the cardinality of the set $\mathcal{T}$.

Note that the model classification listed in

Table 1. Number of variables and constraints in the optimization model (17).

Variables	Type	Number
Active powers	Real	$\left(3n-1\right)p$
Reactive powers	Real	$\left(3n-1\right)p$
Square voltages	Real	$np$
Square currents	Real	$\left(n-1\right)p$
Batteries’ states of charge	Real	$np$
Objective function	Real	2
Total variables	Real	$3\left(3n-1\right)p+2$
Constraints	Type	Number
Active power balance	Equality (affine)	$np$
Reactive power balance	Equality (affine)	$np$
Voltage drops	Equality (affine)	$\left(n-1\right)p$
Powers from Tellegen’s theorem	Inequality conic constraint	$\left(n-1\right)p$
Batteries’ state of charge behavior	Equality (affine)	$np$
Apparent power in batteries	Inequality conic constraint	$np$
Active power cap. in generators and batteries	Inequalities (box-type constraints)	$2np$
Batteries’ state of charge bounds	Inequality (box-type constraint)	$np$
Apparent power in disperse generators	Inequality conic constraint	$np$
Voltage regulation bounds	Inequality (box-type constraint)	$np$
Initial and final state of charges	Equality (affine)	$2n$
Objective functions	Equality (affine)	2
Total constraints	Equalities + Inequalities	$\left(11n-2\right)p+2\left(n+1\right)$

confirms the conic nature of the optimization model given by the power definition at each line and the power outputs in batteries and dispersed sources. In addition, the nature of the decision variables shows that all of those in the studied model are defined in the real domain.

4. Solution Strategy

Due to the fact that this is a multi-objective mathematical model (17) and that convex optimization methods are applicable to single-objective models, in this study, the weighted weights method is used to derive a single objective function [37] in which the weights are varied in order to obtain the Pareto optimal front. The equivalent objective function (i.e., $f_e$) takes the following form:

$$\begin{array}{c}\hfill min{f}_e=\omega \frac{f_1}{F_1^{max}}+\left(1-\omega \right)\frac{f_2}{F_2^{max}},\end{array}$$

(18)

where $\omega$ is the weight factor, associated with each objective function, which ranges from 0 to 1. If $\omega =1$, the global minimum is determined for the objective function associated with greenhouse gas emissions ($F_1^{max}$), and if $\omega =0$, the global optimum will be found for the cost of energy losses during the study period ($F_2^{max}$). The Pareto optimal front is found by starting with $\omega =0$ and increasing its value in 0.05 increments until reaching $\omega =1$ using an iterative process in which the global optimum of the combined problem is found for each value of $\omega$ between 0 and 1. Equation (18) was normalized to combine the two objective functions because they have different units. In addition, in this article, the mathematical model (17) is solved using the CVX optimization package in the MATLAB 2020b programming environment and the SeDuMi and SDPT3 solvers [38].

5. Test System and Characteristics

To evaluate the optimal battery operation method in distribution systems presented in (17), the IEEE 33-node radial distribution system [11] is considered as the test system, which is operated at medium-voltage levels, with a voltage of 12.66 kV at the slack node (bar 1). The structure of this test system is shown in Figure 1. The total active and reactive power consumption at the peak hour of this system is 3715 kW and 2300 kvar, respectively. The data associated with the reactances and resistances of the lines in this test system are available from [11].

The data of the daily demand curves and solar and wind generation profiles are outlined in

Table 2. Solar and wind generation data for each time period.

Per. (h)	PV1 (pu)	PV2 (pu)	WT1 (pu)	WT2 (pu)	Dem. (pu)
1	0	0	0.6303	0.4854	0.4854
2	0	0	0.6194	0.4632	0.4632
3	0	0	0.6098	0.4452	0.4452
4	0	0	0.6050	0.4369	0.4369
5	0	0	0.6122	0.4341	0.4341
6	0	0	0.6411	0.4383	0.4383
7	0	0	0.6927	0.4438	0.4438
8	0	0	0.7395	0.4535	0.4535
9	0	0	0.7779	0.4646	0.4646
10	0	0	0.7887	0.4799	0.4799
11	0	0	0.7671	0.5007	0.5007
12	0	0	0.7479	0.5381	0.5381
13	0	0	0.7287	0.5936	0.5936
14	0	0	0.7371	0.6491	0.6491
15	0	0	0.7731	0.7115	0.7115
16	0.0016	0.0238	0.8031	0.7698	0.7698
17	0.0403	0.1244	0.8271	0.8225	0.8225
18	0.1344	0.2715	0.8523	0.8558	0.8558
19	0.2710	0.4394	0.8788	0.8682	0.8682
20	0.3673	0.6037	0.9064	0.8793	0.8793
21	0.4584	0.7325	0.9328	0.8890	0.8890
22	0.6125	0.8096	0.9520	0.9001	0.9001
23	0.8134	0.8603	0.9640	0.9112	0.9112
24	0.9122	0.9838	0.9700	0.9251	0.9251
25	0.9633	1.0000	0.9748	0.9376	0.9376
26	1.0000	0.9754	0.9784	0.9501	0.9501
27	0.9582	0.7709	0.9832	0.9626	0.9626
28	0.8791	0.8888	0.9880	0.9750	0.9750
29	0.7308	0.6296	0.9940	0.9861	0.9861
30	0.7645	0.5932	0.9988	0.9945	0.9945
31	0.6866	0.4043	1.0000	0.9986	0.9986
32	0.5893	0.3659	0.9964	1.0000	1.0000
33	0.4193	0.2307	0.9892	0.9986	0.9986
34	0.2784	0.2025	0.9652	0.9931	0.9931
35	0.1373	0.1299	0.9244	0.9847	0.9847
36	0.0374	0.0532	0.8607	0.9626	0.9626
37	0.0007	0.0084	0.7743	0.9293	0.9293
38	0	0	0.7251	0.8613	0.8613
39	0	0	0.7167	0.7587	0.7587
40	0	0	0.7167	0.7004	0.7004
41	0	0	0.7251	0.6865	0.6865
42	0	0	0.7263	0.6865	0.6865
43	0	0	0.7179	0.7018	0.7018
44	0	0	0.7095	0.7143	0.7143
45	0	0	0.6987	0.7226	0.7226
46	0	0	0.6915	0.7295	0.7295
47	0	0	0.6867	0.7351	0.7351
48	0	0	0.6831	0.7379	0.7379

and were retrieved from [8]. The renewable generation sources have the following peak values retrieved from [24].

✓

The photovoltaic generator PV₁ is connected to node 13 with a nominal generation of 450 kW, and the photovoltaic generator PV₂ is installed at node 25 with a nominal generation of 1500 kW;

✓

The wind generator WT₁ is installed at node 13 with a nominal power of 825 kW, and the wind generator WT₂ is installed at node 30 with a maximum power generation of 1200 kW.

The IEEE 33-node system relies on three lithium-ion batteries installed, whose characteristics are listed below [39]:

✓

A type-A battery with an energy storage capacity of 1000 kWh and charge/discharge times of 4 h;

✓

A type-B battery with an energy storage capacity of 1500 kWh and charge/discharge times of 4 h;

✓

A type-C battery with an energy storage capacity of 2000 kWh and charge/discharge times of 5 h.

According to [11], the type-A battery is installed at node 14, the type-B battery at node 31, and the type-C battery at node 6, respectively. It is important to note that, in order to preserve the useful life of these batteries, their states of charge must be maintained within a range of 10% to 90%, based on experimental data reported in [40].

To quantify the effects of the variation in demand and alternative energies in energy storage systems, the average rate of CO₂ emitted to the atmosphere by a diesel generation source is set at 612.35 kg/MWh and the energy cost in USD/kWh is set at 0.1390. These data have been set considering an average emission of 1350 lb/MWh and an average energy cost in Colombia for May 2019 of COP/kWh of 479.3389 [8].

6. Computational Validation and Discussion

In this section, the computational validation of the proposed method for the economic–environmental dispatch of energy storage systems in AC distribution networks is presented using a convex multi-objective programming model (see mathematical formulation (17)). The solution of this mathematical model is implemented on a personal computer MD Ryzen 7 3700U (AMD, Santa Clara, CA, USA), 2.3 GHz, 16 GB RAM with 64-bits Windows 10 Home Single Language using the software MATLAB 2020b, CVX convex optimization package, and SeDuMi and SDPT3 solvers.

The CVX package is a mathematical modeling tool compatible with MATLAB software for solving convex optimization problems, making it possible to write objective functions and constraints using the general syntax of MATLAB software. For further details, please refer to [41].

To operate the batteries in the IEEE 33-node test system, two test scenarios are considered: scenario 1 (i.e., E₁), which considers the operation of the batteries with a power factor of 1.0, and scenario 2 (i.e., E₂), which considers the reactive power injection. In a distributed generation, in both scenarios, the batteries operate with a unity power factor. The batteries start and end the daily operating periods charged to 50%, as recommended in [22].

6.1. Scenario 1

Figure 2 shows the Pareto front associated with the BESS operating with a power factor of 1.0.

The results shown in Figure 2 indicate that:

✓

The maximum value of CO₂ gas emissions to the atmosphere, i.e., $F_1^{max}$, is 13.8780 tons/day; whereas the maximum value of the daily cost of energy losses, i.e., $F_2^{max}$, is USD 157.0549;

✓

The minimum values of CO₂ gas emissions to the atmosphere, i.e., $F_1^{min}$, is 6.5502 tons/day; whereas the minimum value of the daily cost of energy losses, i.e., $F_2^{min}$, is USD 132.1042;

✓

The difference between the maximum and minimum values of CO₂ gas emissions is 7.3278 tons/day, which implies a maximum achievable reduction of 52.80%;

✓

The analysis of the extremes of the daily cost of energy losses shows that the difference between its maximum and minimum values is USD 24.9507, which implies a maximum achievable reduction of 15.89%.

Table 3. List of solutions of the Pareto front of Figure 2.

No.	${\mathit{f}}_1$ (tons/day)	${\mathit{f}}_2$ (USD/day)	No.	${\mathit{f}}_1$ (tons/day)	${\mathit{f}}_2$ (USD/day)
1	13.8780	132.1042	7	8.3812	143.4954
2	12.9941	132.3619	8	7.7757	146.8818
3	12.0455	133.2272	9	7.2971	150.0405
4	11.0497	134.8376	10	6.9370	153.0781
5	10.1083	137.1284	11	6.6395	156.0007
6	9.1637	140.2305	12	6.5502	157.0549

presents the values of the objective function for each point represented in the Pareto front of Figure 2.

All solutions presented in Table 3 are optimal because each solution is better than the following in at least one of the two objective functions; that is, the selection of the combination applicable to the distribution system is purely technical and depends solely on the network operator. Nevertheless, the set of solutions 5 to 9 (center of the Pareto front of Figure 2) show an adequate balance between optimizing CO₂ emissions and reducing the daily operating costs. Solution no. 5 reduces emissions by 27.16% and daily operating costs by 12.69%, relative to their maximum values. Solution no. 9 reduces CO₂ emissions by 47.42% and daily operating costs associated with energy losses by 4.47%.

6.2. Scenario 2

In this scenario, the effect of the reactive power injections of the batteries on the Pareto optimal front is evaluated. Figure 3 shows this Pareto front.

The results presented in the Pareto front of Figure 3 show that:

✓

The maximum value of CO₂ gas emissions to the atmosphere, i.e., $F_1^{max}$, is 13.7072 tons/day, whereas the maximum value of the daily cost of energy losses, i.e., $F_2^{max}$, is USD 68.6022. These results imply a reduction of 0.1708 tons/day in greenhouse gas emissions and 85.4527 USD/day in the cost of energy losses;

✓

The minimum value of CO₂ gas emissions to the atmosphere, i.e., $F_1^{min}$, is 6.3317 tons/day, whereas the minimum value of the daily cost of energy losses, i.e., $F_2^{min}$, is USD 44.7601. As in the case of the upper bounds, in this case, daily reductions of 0.2185 tons of CO₂ and USD 87.3441 are achieved;

✓

The main effect of reactive power injection on the multi-objective operation of the BESS is to reduce energy losses in conductors, because this reduction helps to locally offset the consumption of reactants, which translates into currents of lower magnitude in the branches, and, therefore, lower losses.

Table 4. List of solutions of the Pareto front of Figure 3.

No.	${\mathit{f}}_1$ (tons/day)	${\mathit{f}}_2$ (USD/day)	No.	${\mathit{f}}_1$ (tons/day)	${\mathit{f}}_2$ (USD/day)
1	13.7072	44.7601	7	8.0085	56.6499
2	12.7983	45.0283	8	7.4001	59.9527
3	11.8224	45.9312	9	6.9359	63.0869
4	10.7779	47.6243	10	6.5899	65.9802
5	9.8010	49.9735	11	6.3386	68.5226
6	8.8164	53.2112	12	6.3317	68.6023

presents the values of the objective functions for each point represented in the Pareto front of Figure 3.

The comparison of the solutions reported in Table 3 and Table 4 shows that the average reduction in greenhouse gas emissions is 0.2910 tons/day, i.e., over 200 kg of CO₂ per day, and that the reduction of the operational costs of power losses is, on average, USD 87.2598 per day of operation. These results clearly demonstrate the positive effects of reactive power injection through the battery converter, as reported in [8].

6.3. Complementary Analysis

The performance of BESS during one possible operating scenario was illustrated by assuming that the network operator selected solution no. 7, listed in Table 4, whose results associated with the battery variables are shown in Figure 4.

The results displayed in Figure 4 show that:

✓

The batteries are charged during the earlier time periods between 6 a.m. and 10 a.m. (see Figure 4a), periods during which the active power of the battery is negative (see Figure 4b), because this energy is being used to increase the state of charge from 50% to 90%, which corresponds to the upper bound established for this state of charge;

✓

Between 10 a.m. to 6 p.m., the state of charge of the batteries remains stable, with only small variations, which implies that the injected (absorbed) active power is practically nil in this entire operating range;

✓

After 6 p.m., the batteries start discharging the stored energy to the grid, which coincides with a constant decrease in the state of charge until 10 p.m., after which, its state of charge starts to increase again until reaching the final condition of a 50% state of charge, which was set as an operational target. When the state of charge decreases, the active power is positive, which implies a power flow from the battery to the electrical grid (compare Figure 4a,b);

✓

The behavior of the reactive power shown in Figure 4c demonstrates that this variable is independent of the state of charge (see the mathematical model (17)) and that it only depends on the ability of the converter to transfer/absorb power. Accordingly, the battery always delivers reactive power to the grid, and delivers the maximum possible reactive power in the period from 8 a.m. to 11 p.m.

The above variations demonstrate that batteries substantially reduce the cost of energy losses while mitigating greenhouse gas emissions when considering their apparent power injection capabilities because they can perform as dynamic active and reactive power compensators with a range of maneuvers in all four quadrants. This variable performance of the injected reactive power is achieved using control methods (see [21]), enabling better solutions than those achieved when placing a variable pitch capacitor without investing in additional devices.

Figure 5 shows the variation of the states of charge of the three batteries when considering the two extremes of the Pareto front in the case of the unity power factor; that is, the minimization of daily energy losses and greenhouse gas emissions (see points in Figure 2).

The behavior of the states of charge of the batteries located in nodes 6, 14, and 31 in this comparative scenario shows that the main difference between both extreme operating scenarios is identified between periods 16 and 28 for the three batteries because the state of charge is higher when minimizing losses than when minimizing CO₂ emissions. This variation occurs because, during this period, the energy demand markedly increases, which should be partly met by the conventional source. This implies that conventional generation will be lower if the batteries are discharged in this period (batteries injecting active power). Therefore, atmospheric gas emissions will decrease, thus fulfilling the objective of minimizing emissions. In all scenarios analyzed in this study, the batteries always met their operating margins. In other words, the solution will always be feasible.

It is important to mention that the uncertainties of the primary resources, such as wind speed and solar radiation for renewable energy, can be a factor that is essential in the performance of the proposed multi-objective problem performed in this paper. However, using an adequate methodology, this problem can be resolved as presented in [42]. The authors of [42] proposed a mix between a receding horizon control and non-linear autoregressive exogenous model (NARX) to reduce the forecasting errors of the primary resources. This mix consists of using the NARX to predict the wind speed and solar radiation by n-periods and performing the economic dispatch in the $n+1$ period, then recalculating the predictions of wind speed and solar radiation with information of the $n+1$ period.

For implementing the proposed conic model, it is essential to note that two main considerations were made: (i) all of the renewable generation and demand profiles are considered as exogenous inputs in the proposed conic model, which implies that the accuracy between the projected daily operation and the real daily operation results will depend only on the quality of those data; and (ii) the daily battery profiles are sent to the physical location of the batteries using available communication channels, since the proposed optimization methodology corresponds, in other words, to a tertiary controller, which was designed considering a centralized structure.

7. Conclusions and Future Work

In this article, a multi-objective formulation was presented for the problem of optimal battery operation in radial alternating current distribution systems by reformulating a second-order cone programming model to transform the original non-linear problem into a convex equivalent, thereby ensuring that the optimal Pareto front was obtained. For this purpose, the weighted weights method was applied to combine both objective functions so that the multi-objective problem (greenhouse gases vs. cost of energy losses) was solved through an equivalent single-objective programming model. The main advantage of this SOCP formulation is that, through the convex optimization theory, the global optimum is always found using the interior-point or logarithmic barrier method, which only requires transforming the multi-period optimal power flow problem into an equivalent cone programming model by relaxing the hyperbolic relationship between the voltages and the powers in each of the network branches.

Including the dispatch of the reactive power available in the battery’s conversion system, on average, CO₂ emissions are reduced by 200 kg and the system’s operating cost is reduced by more than USD 80 per day of operation, when compared to the unity power factor operating scenario. These results confirmed the benefits of analyzing battery banks as dynamic apparent power compensation banks, which, by taking advantage of electronic power converters, guarantee significant reductions in current flows through the branches, which results in significant environmental and economic benefits, without the need for additional investments.

Based on this article, the following studies can be proposed: (i) to reformulate the proposed SOCP model by including binary variables associated with the optimal location of the batteries, which would generate an MI-SOCP model that can be solved using branching and probing methods combined with interior points, thereby guaranteeing that the global optimum is found; and (ii) to include constraints associated with the useful life of the batteries, energy losses in converters that interface them, and uncertainties of the demand and renewable generation in the optimization model, thus guaranteeing their convex structure for the optimal global solution with the SEDUMI and SDPT3 solvers available for the MATLAB CVX tool.