Dynamic Optimal Power Dispatch in Unbalanced Distribution Networks with Single-Phase Solar PV Units and BESS

1. Introduction

In light of the recent energy crisis, in the midst of the energy transition towards carbon free power systems, the advantages of renewable energy sources (RESs), such as photovoltaic (PV) generators whose growth has been exponential over the last twenty years [1], position them in the center of related policies. However, they suffer from variability and stochasticity of the output power which is their main drawback [2] and a challenge attracting the interest of research community. A promising approach towards a more flexible use of PV sources is their combination with battery energy storage systems (BESSs) to fully exploit the potential of PV units integrated into distribution networks (DNs) [3].

In the last few years, a number of authors have addressed the issues of the optimal power dispatch in DNs with integrated BESSs. These studies have been focused on several main aspects such as: modeling BESSs, defining the BESS function in DNs, i.e. objective functions of the optimization problem, defining control variables, dependent variables and technical constraints, as well as the choice of solution methods.

The BESS is a combination of a voltage source converter and a battery pack, such that the active and reactive power can be controlled independently [4,5,6,7,8,9]. In most research articles, BESS are modeled as a set of equations defining that the stored energy at any given time segment is a function of the stored energy, and charge/discharge power at the previous time segment. Internal losses of BESS are accounted for through the charge/discharge efficiency coefficients or functions. The limits of the reactive power output to be provided by the BESS are also defined in these models.

Depending on the function that the BESS performs in the DN, the objectives for the optimal power dispatch can be summarized in three groups: technical objectives, financial objectives, and multi-objectives. Financial objectives can be related to the maximization of revenues of a RES-BESS based hybrid power plant [10], the reduction in the energy cost from the source grid [11,12], as well as the minimization of distributed generators dispatch costs [13], power generation costs and load shedding costs [5]. Technical objectives are the base objectives and include the static voltage stability improvement [6], the transient voltage stability enhancement [8], the energy loss minimization [9], the voltage fluctuation mitigation [14], and minimization of the voltage deviations at load buses [15]. Researchers often adopt multi-objective functions with different combinations of financial and technical objectives, such as the simultaneous maximization of the total yield from the RES and minimization of the total costs of the energy loss in DNs [4], the simultaneous minimization of the power exchange cost between the DN and transmission network and the penalty cost of the voltage deviation in the DN [7].

In its most general mathematical formulation, the optimal power dispatch in DNs with integrated RES and BESS is a non-linear, non-convex, large-scale, dynamic optimization problem with constraints. There are two approaches to solving it, namely classical optimization methods and metaheuristic methods.

The non-linear optimization models are implemented in the GAMS environment and solved using CONOPT3 solver [4], CPLEX solver [5], and KNITRO solver [11]. The big M method is used in [7] to transform the model into a mixed-integer linear programming problem and solved using a commercial software. The interior-point solver fmincon in MATLAB [9] and the simplex solver through docplex python library [10] are also suggested for the optimal power flow in DNs with the BESS, whereas Giuntoli et al. [13] used the CaADdi framework for interfacing with the IPOPT solver for nonlinear programming. In recent research, some population-based methods, such as PSO, have been proposed to find the optimal power dispatch solution in the DN with the BESS [6,15].

This paper is an expanded and significantly modified version of [12]. The focus of this study is unbalanced DNs with single-phase PV generation and BESS. The main contributions of this study are as follows:

• A dynamic model of the optimal active power dispatch (OAPD) in unbalanced DNs with single-phase PV and BESS is proposed to minimize energy costs from the source network.
• The optimal reactive power dispatch (ORPD) in unbalanced DNs with single-phase PV and BESS is considered. The objective functions are minimization of energy losses and voltage unbalance at three-phase load buses, and control variables are reactive powers of solar PV inverters and BESS.
• A scenario based approach encompassing the MSC method and the simultaneous backward reduction technique is proposed for modeling the uncertainty of PV generation and load in the DN.
• Application of a novel hybrid PSOS-CGSA algorithm to solve the optimal active-reactive power dispatch problems is presented.
• Defining a modified IEEE 13-bus feeder for evaluating the solvability and applicability of the proposed approach is suggested.

The rest of the paper is organized as follows. In Section 2, a scenario-based model is formulated for the optimal dispatching of active power in DN by operating mode control of BESSs. The problem of optimal reactive power dispatch in unbalanced DNs with single-phase PV and BESS inverters is elaborated in Section 3. The solution method based on the hybrid PSOS-CGSA metaheuristic algorithm is briefly presented in Section 4. The simulation results are discussed in Section 5, whereas the main conclusions are listed in Section 6.

2. Active Power Dispatch by BESS

Battery operation mode is determined by two basic parameters. These are charge/discharge power and state of charge (SoC). Since the battery capacity is limited, the SoC of the battery must be viewed as a dynamic quantity, as follows [10]:

$$So{C}_{\varphi }\left(t\right)=So{C}_{\varphi }\left(t-\Delta t\right)+\Delta So{C}_{\varphi }$$

(1)

(2)

where $C_{bat,\varphi }$ is the total capacity of the battery; $P_{bat,ch}$ and $P_{bat,dch}$ are charging and discharging power of the battery, respectively; ${\eta }_{ch}$ and ${\eta }_{dch}$ are charging and discharging efficiency, respectively; Δt is the time segment (for example 1 h), and ϕ denotes the phase, ϕ$\in$(a,b,c).

Battery power and SoC must always be within defined limits, as follows:

$$0\le P_{bat,ch,\varphi }^{}\left(t\right)\le P_{bat,ch\max }^{}$$

(3)

$$0\le P_{bat,dch,\varphi }^{}\left(t\right)\le P_{bat,dch\max }^{}$$

(4)

$$So{C}_{\min }^{}\le So{C}_{\varphi }\left(t\right)\le So{C}_{\max }^{}$$

(5)

$$So{C}_{\varphi }\left(t=0\right)=So{C}_{\varphi }\left(t=T\right)$$

(6)

$$P_{bat,ch,\varphi }^{}\left(t\right)\cdot P_{bat,dch,\varphi }^{}\left(t\right)=0$$

(7)

where $P_{bat,ch\max }^{}$ and $P_{bat,dch\max }^{}$ are maximum charging and discharging powers of the battery, respectively; $So{C}_{\min }^{}$ and $So{C}_{\max }^{}$ are the predefined minimum and maximum charge levels, respectively; T is a time horizon, e.g. 24 h.

Battery charge/discharge power ($P_{bat}$) and SoC are mutually dependent variables. Taking into account the variability of the electricity price, their relationship can be expressed as follows [12]:

(8)

where $C_g^{}(t)$ is the day-ahead electricity price from the source grid at time t, $C_{gav}^{}$ is the average electricity price during the day; a_ϕ(t) is a coefficient within the range [0, 1].

The coefficient a_ϕ(t), which defines the power of the BESS, can be determined in different ways. In [12], the coefficient a_ϕ(t) is calculated based on the relation of the electricity price at the time t and the mean electricity price during a day. A similar approach was adopted in [16], where the coefficient a_ϕ(t) is calculated as a function of the difference in the electricity price from the predefined lower price limit for discharging and upper price limit for charging. However, the charging/discharging powers of the BESS determined in this way do not guarantee the optimal dispatch of the active power according to the adopted objective function. Therefore, in this paper, the coefficient a_ϕ(t) is considered as a control variable, whose value optimizes to ensure the minimization of a given objective function.

The basic equation of the active power dispatch problem is the power balance constraint (9), which must be satisfied at each time segment t:

$${\displaystyle \sum _{\varphi \in \left(a,b,c\right)}{P}_{g,\varphi }^{}\left(t\right)}={\displaystyle \sum _{\varphi \in \left(a,b,c\right)}\left[P_{L,\varphi }^{}(t)-P_{PV,\varphi }^{}\left(t\right)+P_{bat,\varphi }^{}\left(t\right)\right]}$$

(9)

where P_g,ϕ(t) is the active power from/to the source grid at time t in phase ϕ; P_L,ϕ(t) is the total active power of the load at time t in phase ϕ; P_PV,ϕ(t) is the total active power generation of PV units at time t in phase ϕ; P_bat,ϕ(t) is the active power of the battery at time t in phase ϕ; ϕ$\in$(a,b,c).

2.1. Objective Function for OAPD

By controlling the operating mode of the battery, i.e. the charging/discharging power and SoC, the active power flows in the distribution network can be significantly influenced. In this way, the exchange power with the source power grid is directly affected, as well as other parameters of the distribution networks, such as power losses, voltage profile, etc. The operating costs of the distribution network mostly depend on the value of the energy from/to the source grid. Considering various PV generation and load demand scenarios with their probabilities, the objective function for OAPD can be formulated as the minimization of the expected cost for electricity, as follows:

(10)

(11)

(12)

where ω and Ω are the current scenario and the total number of scenarios, respectively; p(ω) is the probability of scenario ω; $C_g^{\omega }$ and $C_{bat}^{\omega }$ are the total cost for electricity from the source grid and operating cost of the battery in scenario ω, respectively; c_g and$c_{bat}$ are the electricity price from the source grid and the unit operating cost of the battery, and T is the considered period (i.e. 24 h).

Two assumptions should be kept in mind here, namely:

• The retail electricity prices have the same daily profile as the wholesale electricity prices and fluctuate hourly.
• The DN transacts with the wholesale markets through some energy agent such as an aggregator.

The power output of a PV source is determined by the solar irradiation at a given moment, and its value cannot be influenced. Therefore, the cost of the PV generation is not taken into account in the objective function (10), because it cannot be optimized with controlling the charge/discharge power of the battery.

For stable operation, constraints (3-7) and (9) must be satisfied. The output power of the PV generator is limited by the maximum output power, i.e. the nominal power. Also, P_g,ϕ(t) must be restricted by minimum and maximum values.

2.2. Modeling of Uncertainties

2.2.1. PV Generation Modeling

The power of the PV source as a function of solar irradiance can be expressed as follows [16]:

$$P_{PV}\left(S\right)=\{\begin{array}{cc}{P}_{pvn}\frac{S^2}{S_{stc}{R}_c}& \mathrm{for}S<R_c\\ P_{pvn}\frac{S}{S_{stc}}& \mathrm{for}S\ge R_c\end{array}$$

(13)

where P_pvn is the nominal power of the PV source; S is the solar irradiance on the PV module surface; S_stc is the solar irradiance under standard test conditions; R_c is a certain irradiance point.

To model the stochastic nature of the solar irradiance, it is most convenient to use the Beta probability density function (PDF):

$$f_b\left(S\right)=\{\begin{array}{c}\begin{array}{l}\frac{\Gamma \left(\alpha +\beta \right)}{\Gamma \left(\alpha \right)\Gamma \left(\beta \right)}{S}^{\left(\alpha +1\right)}{\left(1-S\right)}^{\left(\beta -1\right)},\\ \mathrm{for}0\le S\le 1,\alpha \ge 0,\beta \ge 0\end{array}\\ 0,\mathrm{otherwise}\end{array}$$

(14)

where f_b(S) is the Beta function of the solar irradiance S; α, β are the shape parameters of the Beta distribution; and Γ is the Gamma function.

Based on the data of the long-term measurement of the solar irradiance in a given area, characteristic daily diagrams of the solar radiation can be estimated, with statistical indicators, i.e. mean value (µ_S) and standard deviation of the solar irradiance (σ_S) in a given time interval t (e.g. an hour). Using mean value and standard deviation of the solar irradiance, the shape parameters of the Beta PDF can be calculated for a time interval t, as follows:

(15)

$$\alpha =\frac{{\mu }_S\cdot \beta }{1-{\mu }_S}$$

(16)

2.2.2. Load Modeling

The normal (Gaussian) probability function is most commonly used to model load uncertainty in distribution systems. Generally, the load level is assumed to be a random variable (L) following the same PDF within each hour of a given daily load diagram [2].

$$f_n\left(L\right)=\frac{1}{\sqrt{2\pi }{\sigma }_L}{e}^{-\frac{{\left(L-{\mu }_L\right)}^2}{2{\sigma }_L^2}}$$

(17)

$$L={\mu }_L+\sqrt{2}{\sigma }_L{\mathrm{erf}}^{-1}\left(2r-1\right)$$

(18)

where μ_L and σ_L are the mean value and standard deviation of L, respectively; r is a random number in [0, 1]; and erf(·) and erf^-1(·) are the error function and the inverse error function, respectively.

2.2.3. Scenario Generation

The Monte Carlo simulation (MCS) method is suitable to modeling the uncertainties of solar irradiance and load level. Starting from probability density functions for solar irradiance and load, whose parameters are determined based on historical data, eg previous years, the MSC method can generate a large number of scenarios for day-ahead profiles of solar irradiance and load levels. Although a larger number of scenarios (e.g. several thousand) provides a more faithful modeling of the uncertainty of the output power of the PV source and the load level, on the other hand, it means a higher computation burden. In distribution network operation control applications, such as optimal power flow, methods with high accuracy and low computational burden are needed.

As a compromise between a good approximation (modeling) of the uncertainty of stochastic input variables (solar irradiance and load) on the one hand, and a significant reduction in calculation time on the other hand, a scenario reduction strategy aggregating similar scenarios can be applied. In doing so, it must be ensured that the metric of probability distribution of the reduced set must be close enough to the original set of scenarios. There are several algorithms for reducing the number of scenarios, such as k-means clustering, backward reduction, and fast forward selection. In this work, the simultaneous backward reduction algorithm [17,18,19] was employed to reduce a ternary scenario tree for the daily PV generations and load profiles.

3. Optimal Reactive Power Dispatch

In this paper, it is assumed that PV arrays and batteries are connected to unbalanced DNs via single-phase converters. Depending on the nominal apparent power of the converter and the active power level of the corresponding PV source or battery, they have the ability to control the generation/absorption of reactive power, that is, they can be used as an additional regulation resource in DNs [12,20], as illustrated in Figure 1.

Therefore, for DNs with PV units and BESS, the problem of ORPD can be posed as a problem of determining the optimal reference values of the available reactive powers of the inverters and the slack bus voltage magnitudes with certain objective function, under specified active power outputs of all PV and BESS inverters representing the solution of the active power dispatch as explained in Section 2.

3.1. Objective Function for ORPD

3.1.1. Active Power Loss Minimization

The level of active power losses is one of the most important indicators of the efficiency of the distribution network management, because they represent the direct costs for a distribution company. Therefore, the minimization of active power losses is one of the most common and important objective functions in OPF studies used by distribution system operators. Considering different (reduced) scenarios for PV generation and load profiles with their respective probabilities, the objective function for ORPD can be formulated as a minimization of the expected daily energy loss, as follows:

$$Fobj1=\min \left[{\displaystyle \sum _{\omega =1}^{\Omega }{W}_{loss}^{\omega }p\left(\omega \right)}\right]$$

(19)

$$W_{loss}^{\omega }={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{\begin{array}{l}i=1\\ \varphi \in a,b,c\end{array}}^{NB}{R}_{i,\varphi }{\left(I_{i,\varphi }^{\omega }\left(t\right)\right)}^2}}$$

(20)

where $W_{loss}^{\omega }$ is the total energy loss in scenario ω; NB is the number of branches in the DN; $R_{i,\varphi }$ is the resistance of the branch i in phase ϕ; $I_{i,\varphi }^{\omega }$ is the current magnitude through the ith branch in phase ϕ for scenario ω.

3.1.2. Voltage Unbalance Minimization

High penetration of single-phase PV systems and relative variations in per-phase loading lead to an increase in voltage imbalance in distribution networks. It is well known that significant voltage unbalance may lead to higher network losses, increased heating of network components such as transformers, as well as costly damage or derating of three-phase motors [21]. In this paper, the phase voltage unbalance rate (PVUR) is used for the voltage unbalance metric. The PVUR is calculated at each three-phase bus using the line-to-ground voltage magnitudes V_a, V_b, V_c [21]:

$$\begin{array}{l}PVUR(\%)=\frac{\Delta V_P^{\max }}{V_P^{avg}}\cdot 100,\\ \mathrm{where}{V}_P^{avg}=\frac{V_a+V_b+V_c}{3},\\ \Delta V_P^{\max }=\max \left\{\left|V_a-V_P^{avg}\right|,\left|V_b-V_P^{avg}\right|,\left|V_c-V_P^{avg}\right|\right\}\end{array}$$

(21)

According to the IEEE standard 141-1993, the maximum value of PVUR must not exceed 2% for distribution networks. To ensure that the PVUR value at any three-phase bus does not exceed this limit, the minimization of the expected maximum of PVUR in the DN can be defined as an objective function:

$$Fobj2=\min \left[{\displaystyle \sum _{\omega =1}^{\Omega }PVU{R}_{\max }^{\omega }p\left(\omega \right)}\right]$$

(22)

where $PVU{R}_{\max }^{\omega }$is the maximum value of the PVUR in the DN for scenario ω.

3.2. Constraints

3.2.1. Equality Constraints

The equality constraints for the ORPD problem are typical nonlinear power flow equations for three-phase unbalanced DNs [22,23].

3.2.2. Inequality Constraints

There are two types of inequality constraints: functional operational constraints such as slack bus active power limits (23), load bus voltage limits (24), and branch flow limits (25); and feasibility region constraints defined by inverter reactive power limits (26), and slack bus voltage limits (27).

$$P_{grid,\min }^{\varphi }\le P_{grid}^{\varphi }\left(t\right)\le P_{grid,\max }^{\varphi }$$

(23)

$$V_{i,\min }^{\varphi }\le V_i^{\varphi }\left(t\right)\le V_{i,\max }^{\varphi },i=1,\dots ,NL$$

(24)

$$S_i^{\varphi }\left(t\right)\le S_{i,\max }^{\varphi },i=1,\dots ,NB$$

(25)

$$-\sqrt{{\left(S_{in,i}^{\varphi }\right)}^2-{\left(P_{in,i}^{\varphi }\left(t\right)\right)}^2}\le Q_{in,i}^{\varphi }\left(t\right)\le \sqrt{{\left(S_{in,i}^{\varphi }\right)}^2-{\left(P_{in,i}^{\varphi }\left(t\right)\right)}^2}$$

(26)

$$V_{0,\min }^{\varphi }\le V_0^{\varphi }\left(t\right)\le V_{0,\max }^{\varphi }$$

(27)

where $P_{grid}^{\varphi }\left(t\right)$, $P_{grid,\min }^{\varphi }$ and $P_{grid,\max }^{\varphi }$ are the active power from/to the source grid at time t in phase ϕ, minimum and maximum power from/to the source grid, respectively; $V_i^{\varphi }\left(t\right)$, $V_{i,\min }^{\varphi }$, $V_{i,\max }^{\varphi }$ are the voltage magnitude at the load bus i in phase ϕ, minimum and maximum limits of the voltage magnitude at load buses, respectively; $S_i^{\varphi }\left(t\right)$ and $S_{i,\max }^{\varphi }$ are the apparent power value and maximum apparent power value of the ith branch in phase ϕ, respectively; $Q_{in,i}^{\varphi }$ is the reactive power injection of the inverter connected to phase ϕ at the ith bus, $S_{in,i}^{\varphi }$ and $P_{in,i}^{\varphi }\left(t\right)$are the maximum (rated) apparent power and the active power at time t of the PV (or BESS) inverter connected to phase ϕ at the ith bus, respectively; $V_0^{\varphi }\left(t\right)$, $V_{0,\min }^{\varphi }$ and $V_{0,\max }^{\varphi }$ are the voltage magnitude value, minimum and mаximum limits of the voltage magnitude at the slack bus in phase ϕ, respectively.

4. Solution Method

The hybrid PSOS-CGSA algorithm introduced in [23] was applied to effectively solve the above defined optimal active/reactive power dispatch. This is a metaheuristic population-based optimization algorithm constructed by the hybridization of particle swarm optimization with sigmoid-based acceleration coefficients (PSOS) [24] and chaotic gravitational search algorithm (CGSA) [25]. Its performance in terms of simplicity, efficiency and robustness has been demonstrated in solving various complex optimization problems such as distribution system state estimation [23], dynamic economic dispatch [26], and optimal power flow in DNs [12].

4.1. Application of PSOS-CGSA for Active Power Dispatch

In the optimal active power dispatch problem presented in Section II, the objective function is defined by equation (10). The control variables are coefficients a_ϕ(t) which define the power of the BESS according to (8). Therefore, a search agent can be written as:

$$x_i\left(t\right)=\left[a_{a,i}\left(t\right),a_{b,i}\left(t\right),a_{c,i}\left(t\right)\right]i=1,\mathrm{..},N$$

(28)

where N is the total number of search agents (population size), a_ϕ,i(t) is the value of control variable a within search agent i for phase ϕ at iteration t.

A summary of the proposed algorithm for solving the active power dispatch is represented by the flowchart in Figure 2.

4.2. Application of PSOS-CGSA for Reactive Power Dispatch

The vector of control variables in the problem of ORPD consists of the reactive powers of the PV and BESS inverters and the slack bus voltage magnitudes, and can be expressed as follows:

$$x_i\left(t\right)=\left[Q_{i,1}^{\varphi }\left(t\right),\mathrm{...}{Q}_{i,j}^{\varphi }\left(t\right),\mathrm{...},Q_{i,N_i}^{\varphi }\left(t\right),V_0^{\varphi }\left(t\right)\right]i=1,\mathrm{..},N$$

(29)

where $Q_{i,j}^{\varphi }$ is the reactive power of the inverter at the jth bus in phase ϕ, N_i is the total number of buses in which inverters are connected, $V_0^{\varphi }$is the voltage magnitude at the root bus of phase ϕ$\in$(a,b,c).

The application of PSOS-CGSA in solving the problem of optimal reactive power dispatch is presented with the flowchart in Figure 3.

5. Simulation Results

The proposed approach for the day-ahead optimal dispatching of active and reactive powers was simulated on a modified IEEE 13-bus feeder [27]. In the modified IEEE 13-bus feeder, the single-phase PV sources are connected to buses 634, 645, 646, 611, 675, 692 and 652, as shown in Figure 4. Single-phase PV sources, each with the same rated power of 250 kW, are connected to the grid via single-phase inverters. Single-phase BESSs, each with a rated power of 0.9 MW and a capacity of 3.6 MWh are connected at bus 680 in all three phases via single-phase converters rated apparent power 1.1 of the rated power of the batteries. Data on PV and BESS parameters are given in Table 1 and Table 2, respectively.

The daily load profile is assumed to be with the hourly mean values tabulated in Table 3 with the same standard deviation of 10% of the corresponding mean values. The simulation was performed for the daily solar irradiance profile with the hourly mean value and standard deviation as given in Table 4. The day-ahead electricity price (EP) data are adopted as in Table 5. The maximum value is estimated according to the electricity price for Serbia on the spot market in January-February, 2023.

5.1. OAPD Results

As explained in Section 2, using MCS, there are 2000 solar irradiace and load scenarios generated, as shown in Figure 5. By applying the simultaneous backward reduction technique, these 2000 scenarios were reduced to 10 scenarios shown in Figure 6.

By applying the proposed approach for the optimal active power dispatch, the expected total operating costs of each scenario were determined, and the results are shown in Table 6.

The optimal value of the expected daily cost for electricity is 453.16 (€). Acording to (10), this value is obtained as sum of the costs in each reduced scenario multiplied by the corresponding probability. In view of this, the expected day-ahead charging/discharging schedule of the BESS was obtained as shown in Figure 7 and Figure 8.

To fair comparison, Figure 9 shows the daily costs for electricity in the following cases: (i) without the BESS integrated into the DN - base case, (ii) with the BESS in the DN applying the charge/discharge schedule approach presented in [12], and (iii) with the BESS in the DN using the proposed OAPD algorithm. The results in Figure 9 clearly indicate that the proposed PSOS-CGSA-based approach for OAPD leads to a large reduction in costs for electricity compared to the base case and the approach for OAPD in [12].

5.2. ORPD Results

Here, the optimal values of 21 control variables should be determined, as follows:

$$x=\left[\begin{array}{l}{Q}_{634}^a,Q_{634}^b,Q_{634}^c,Q_{645}^b,Q_{645}^c,Q_{646}^b,Q_{646}^c,Q_{692}^a,Q_{692}^b,Q_{692}^c,\mathrm{...}\\ Q_{675}^a,Q_{675}^b,Q_{675}^c,Q_{611}^c,Q_{652}^a,Q_{680}^a,Q_{680}^b,Q_{680}^c,V_{650}^a,V_{650}^b,V_{650}^c\end{array}\right]$$

The simulation studies performed for the following two test cases:

Case 1: Minimization of the active power losses (Fobj1)

Case 2: Minimization of the phase voltage unbalance rate (Fobj2)

The voltages at the load buses must be within [0.94 -1.06] (p.u.). The voltage of the root bus (V₀) can be changed in the range [0.9 -1.1] (p.u.), while the limits of the reactive power of the converter are variable, according to (26).

The results are given in Table 7 and Table 8. In the base case, the PV sources and BESS operate with unity power factor, and the root bus voltage is equal to 1 p.u. Table 7 and Table 8, respectively, show the expected values of total active energy losses (W_loss) and maximum values of phase voltage unbalance rate (PVUR_max) for all ten reduced scenarios.

Applying the proposed procedure for ORPD achieves a reduction of expected energy losses by 22% compared to the base case, while the expected value of the maximum PVUR drops from 9.1% to 1.91%. There are significant energy savings, and it is ensured that the PVUR index level is below the maximum limit of 2% for all three-phase buses in the DN.

Based on the results shown in Table 7 and Table 8, it can be seen that in the case of the minimization of active energy losses (Fobj1) the best scenario is 9, and in the case of the minimization of the voltage unbalance (Fobj2) the best scenario is scenario 1. Figure 10 and Figure 11 show the optimal values of the control variables for these best scenarios.

The voltage profiles at the bus 671 for the best scenario in Case 1 and the best scenario in Case 2 are shown in Figure 12 and Figure 13, respectively.

In Case 1, ORPD reduces the active energy losses compared to the base case, while the voltages are within the permissible limits, as shown in Figure 12. The phase voltage magnitudes are much closer to each other in Case 2 (Figure 13) than in Case 1 (Figure 12). Clearly, this is as a consequence of the PVUR minimization in Case 2.

5.3. Evaluation of PSOS-CGSA

Тhe efficiency of the proposed solution algorithm (PSOS-CGSA) was confirmed by a comparison with well-established metaheuristic algorithms PSO, GSA and theirs recent improved versions PSOS [24] and CGSA [25]. Algorithm parameters displayed in Table IX, are adopted based on their default values in the cited references. The population size and maximum iteration number are set as: N=50 and it_max=200 for all case studies. All algorithms were employed to solve the problem of OAPD under the same terms and input data.

The daily costs for electricity obtained for all 10 reduced scenarios are shown in Figure 14. It is clear that the proposed hybrid PSOS-CGSA method finds better solutions for the OAPD problem, i.e., provides lower costs in each of 10 scenarios compared to PSO, GSA, PSOS and CGSA algorithms.

Table 9. Algorithms parameters.

PSO	C1=2; C2=2; wmax=0.9; wmin=0.4
GSA	G0=1; α=20; p=2; K0=1
PSOS	C1i=0.5; C1f=2.5; λ=0.0001; wmax=0.9; wmin=0.4
CGSA	G0=1; α=20; wmax=0.1; wmin=1e-10; Sinusoidal map
PSOS-CGSA	C1i=0.5; C1f=2.5; λ=0.0001; G0=1; α=20; Sinusoidal map

Table 9. Algorithms parameters.

PSO	C1=2; C2=2; wmax=0.9; wmin=0.4
GSA	G0=1; α=20; p=2; K0=1
PSOS	C1i=0.5; C1f=2.5; λ=0.0001; wmax=0.9; wmin=0.4
CGSA	G0=1; α=20; wmax=0.1; wmin=1e-10; Sinusoidal map
PSOS-CGSA	C1i=0.5; C1f=2.5; λ=0.0001; G0=1; α=20; Sinusoidal map

6. Conclusion

In this paper, a two-stage approach for the dynamic day-ahead optimal active-reactive power dispatch in unbalanced DNs with high penetration of single-phase solar PV systems and BESS has been proposed. The simulation results conducted on a modified IEEE 13-bus test system indicate the conclusions that can be summarized as follows:

• The proposed OAPD ensures minimum costs for the electricity from the source grid for all scenarios generated by the MSC method and the simultaneous backward reduction algorithm.
• It has been shown that single-phase PV inverters and BESS can serve as additional reactive power management resources.
• The proposed ORPD approach enables the reduction of active power losses in the DN and reduction of the voltage unbalance on the three-phase load buses below the permitted limit.
• The proposed PSOS-CGSA solution technique provides better solutions for the optimal power dispatch problems compared to other well-established metaheuristic algorithms such as PSO, GSA, and modified PSOS and CGSA.