An Optimal Power Flow Algorithm for the Simulation of Energy Storage Systems in Unbalanced Three-Phase Distribution Grids

An optimal power flow algorithm for unbalanced three-phase distribution grids is presented in this paper as a new tool for grid planning on low voltage level. As additional equipment like electric vehicles, heat pumps or solar power systems cause unbalanced power flows, existing algorithms have to be adapted. The focus in this paper is on the simulation of energy storage systems in unbalanced systems. The used grid model, constraints, objective function and solver are explained in detail. A validation of the algorithm using a commercial power flow tool is done. Additionally, two exemplary optimizations are performed to show possible applications for this tool.


I. INTRODUCTION
The planned reduction of carbon dioxide emissions leads to a transition of the energy system towards a renewable energy based generation in many countries worldwide. Several renewable generation units like solar power systems are connected to the low voltage level. Additionally, new loads such as electric vehicles or heat pumps are integrated into the power grid. These systems are mainly connected to the low voltage level as well. The new equipment causes wider fluctuations in the power flow and hence the planning processes are getting more complicated for low voltage grids.
Optimal power flow (OPF) algorithms are known as a planning tool on transmission level. Taking into account the developments mentioned above, applications for this algorithm on distribution level occur. The main difference is that on a distribution level, the assumption of balanced power flows is not given in general. Adapted algorithms have already been presented in literature such as three-phase optimal power flow (TOPF) [1] [2] [3] and used for applications as Conservation Voltage Reduction [4] or Voltage Unbalance Mitigation [5].
In this paper, an algorithm is developed that can be used to simulate energy storage systems in low voltage grids. The focus is hereby on a grid-friendly operation of the storage system [6] in the provided examples. Therefore, a dynamic optimal power flow algorithm like in [1] [7] is necessary. In comparison to other TOPF like [1] [4] [8], that are using a sequence of algorithms internally, the algorithm consists only of a single optimization problem that is build up and solved. As explained later in the paper, the optimization problem size increases in comparison to an algorithm for balanced grids significantly. Nevertheless, simplification techniques that are used in TOPF formulations in literature (e.g. linearization in [9] [10] or relaxation techniques in [11] [12]) are not used as the problem can be solved using the presented solver without approximations. Until now, only few TOPF formulation are also dealing with energy storage systems (e.g. [13]). The novelty of the algorithm presented here is that energy storage systems are considered without the use of simplification techniques.
In Chapter II, the approach for modelling the low voltage grid is presented. In Chapter III and IV, the basic formulation of the optimization problem is explained as well as the used solver. In Chapter V, VI and VII, the exact cost function, optimization variables and constraints are introduced. In Chapter VIII, a validation using a commercial power flow software is done. In Chapter IX, the results of an exemplary optimization including a battery storage system are explained.

II. GRID MODEL
The basic approach to model a grid in most existing OPF algorithms can be seen in Figure 1 for a grid consisting only of one generator, one line and one load. The generator is inserting power to the grid that is flowing over a line to the load, where the energy is consumed. The algorithm presented in this paper is based on [6] and [7], where this approach was used. In Germany, the energy system is designed as a three-phase system. Assuming a balanced power flow on all three phases, it is sufficient to calculate the power flow on only one of the three phases. As the algorithm described in this paper is adjusted especially for the modelling of distribution grids, the assumption of a balanced power flow is not made here in comparison to [6] and [7]. In a typical low voltage feeder in Germany, between ten and 50 households are connected. This number is not high enough to assume a balanced power flow for statistical reasons. Therefore, all three phases have to be modelled independently. Furthermore, the power flow on the neutral line has to be taken into account as an additional voltage drop arises there. The vectorial sum of the current flowing through the different phases of a load is flowing back on the neutral line. Hence, it is necessary to model the load as an impedance where a current is flowing through.
In comparison to the single phase OPF, there are in total four slack buses. One for each phase and on the neutral line an additional generator has to be installed which acts as a slack for the neutral line. This concept is shown in Figure 2. The grid model is used to build up the nodal admittance matrix y.
For the impedances of the lines (Z L1 , Z L2 , Z L3 ,and Z N in Figure 2), the positive-sequence impedances of the lines are used. To simplify the grid model, the negative-and zero-sequence system are not taken into account. Hence, the coupling between the different phases is only partly included. As impedances are used to represent a load, the input data has to be processed to calculate the value of the impedance. In literature, there exist different models for loads. In most OPF calculations, constant power loads are used.
As the power is fixed for a constant power load, the impedance only depends on the voltage as shown in Equation 1. As the voltage is an optimized variable, the impedance of the load changes in each iteration of the optimization. As these impedances are part of further constraints (see Equation  37 and 38), this leads to a situation where the optimization problem is not well-defined using constant power loads. For the test cases shown in this paper, the algorithm is converging assuming constant power loads (e.g. in Chapter VIII), but more iterations are needed.
If the impedance has a constant value, the optimization problem can be formulated correctly. Hence, constant impedance loads are assumed. In this approach, it is assumed that the nominal power is consumed at nominal voltage. So the impedance is calculated as follows III. DYNAMIC OPTIMAL POWER FLOW ALGORITHM Optimal power flow (OPF) problems can be described in the following form: F(x) is the so-called cost function, which is minimized while the equality constraints g(x) as well as the inequality constrains h(x) have to be fulfilled. A dynamic optimal power flow algorithm with a horizon of T is used in this paper to be able to include equipment, whose current condition is dependent on states in previous or future time steps, e.g. any kind of storage system. Therefore the state vector x contains optimization variables for all time steps.
Additionally, the cost function sums up to IV. SOLVER The Primal-Dual Interior Point Method (PDIPM) is used to solve this problem. The Langrangian multipliers µ (inequality constraints) and λ (equality constraints) are assigned. The slack variables Z are introduced to transform the inequality constraints to equality constraints. The slack variable Z is weighted with the barrier coefficient γ to keep the inequality constraints away from zero in the first iterations. Finally, the Langrangian L is built up as follows: The first order optimality conditions for the optimization problem are satisfied when the partial derivatives of the Lagrangian above are all set to zero. In the PDIPM, the first order optimality conditions are solved using Newton's method. Therefore, first and second derivative of the constraints as well as the cost function have to be determined. A detailed explanation of PDIPM can be found in [14].

V. OPTIMIZATION VARIABLES
All optimization variables for one time step are saved in the vector x. The number of optimization variables are determined in the following section for an exemplary grid containing only two points (see Figure 2) for a single time step. At first, vector x contains the real (e t X,Y ) and imaginary (f t X,Y ) parts of the voltages on Phase Y at all Points X. For the exemplary grid, these are 16 variables (2 parts per voltage · 2 points · 4 voltages per point). Additionally, the active power P t G,N as well as the reactive power Q t G,N of all generators are part of the vector x. N is a number assigned to each generator. As there are four generators, there are in total 8 optimization variables. To include a battery storage system number N, the energy capacity E t S,N of the system is one variable in the vector x. Besides, the active charging P t SC,N,Y , active discharging power P t SD,N,Y and reactive power Q t S,N,Y per Phase Y needs to be considered. Including one battery system leads to 10 optimization variables.
In total, optimization vector x contains 24 variables without and 34 variables including the battery system. In comparison, the optimization vector x would contain 6 variables without and 10 variables including the battery system for a traditional balanced power flow formulation.

VI. COST FUNCTION
For the cost function, costs are taken into account using cost factors c. These factors are multiplied with a corresponding power and the chosen temporal resolution. In Equation 9, the cost function is shown where the energy of the generators is priced. These generators represent the power consumption from the local utility. For all cost factors, a price of 28 ct/kWh, which is representative in Germany, is used in this paper. The optimization algorithm consequently optimizes the cost for energy in the regarded grid. (9) VII. CONSTRAINTS

A. Linear Equality Constraints
1) Voltage at the Slack Bus: The complex voltage is divided in two parts, the real and imaginary part.
The magnitude of all voltages at the slack bus is set to 1 p.u. using linear equality constraints. In a three-phase system, there is a 120 degree phase shift between the phases. This phase shift is fixed at the slack bus here to model the feed-in from the medium voltage grid.
Additionally, the magnitude and phase of the voltage at the slack bus on the neutral line is set to zero.
B. Linear Inequality Constraints 1) Power Limits of the Generator: In general, for the active and reactive power of each generator an upper and lower power limit can be defined.
In the test grid, there are only generators at the slack bus. In a low voltage grid, the slack bus represents the connection to the MV/LV transformer. Hence, the limits are set to the power limit of the MV/LV transformer in the modelled grid. In Equation 23 and 24, this is shown exemplary for Phase 1 and hence Generator 1. As the voltage on the neutral phase is set to zero (see Equation 17 and 18), the power of the generator on Phase N does not have to be constrained as no power is flowing there.
2) Energy Capacity Limits of the Storage: The energy capacity of the storage system is limited between an empty storage and the maximum energy capacity E S,N,max of the storage system.
3) Power Limits of the Storage: Due to the construction of a storage system, also the charging and discharging power of the storage is limited.

C. Nonlinear Equality Constraints 1) Storage Equation:
The energy capacity in the current time step E t S,N depends on the energy capacity of the previous time step E t-1 S,N and the energy difference caused in the current time step ∆E t S,N .
The charging efficiency η c and discharging efficiency η d are taken into account as a constant. Alteration of battery parameters through degradation over the battery lifetime are not taken into account.

2) Power Flow Equation:
At each point in the grid, the power flowing into the point from lines as well as the power generated or consumed at these points have to be balanced. This is taken into account using the following constraints Bus is thereby a vector containing the sum of the active power generated or consumed at all point X. P t Line is a vector containing the sum of the active power flowing from a Point X to other points on a line and calculated as follows in Equation 37, where • denotes point-wise multiplication. Q t Bus and Q t Line are defined accordingly.
The nodal admittance matrix y can be separated into a real and an imaginary part as denoted in Equation 39. Equation 37 and 38 are determined using Equation 40. The vector v t contains all voltages v t X,Y for time step t.
D. Nonlinear Inequality Constraints 1) Power Limit of the Lines: The power flow on the branches S t Br is limited using a maximum apparent power S t Br,Max that is allowed to flow over each branch. In the algorithm, the power is limited at each end of the lines.
2) Voltage Limits: For the voltage at each point, upper as well as lower voltage limits can be defined. For the three-phase algorithm, it is important to mention that the valid voltage at consumer level is the difference between phase voltage and neutral line voltage. Hence, for each point on Phase 1, 2 or 3, the difference between phase and neutral line voltage is constrained between a minimum allowed voltage v Min and a maximum allowed voltage v Max .

E. Total Number of Constraints
The total number of constraints for a balanced and unbalanced formulation with and without storage system of the grid in Figure 2 can be seen in Table I. This algorithm is designed as planning tool. Hence, real-time capability is not necessary, but nevertheless a fast computation time is advantageous. Using a standard state of the art computer, the calculations presented in Chapter VIII and IX can be performed in less than one minute. The algorithm is validated using the commercially available software environment DIgSILENT PowerFactory 2019. This is a load flow software, so an optimization of equipment is not possible and therefore this application can not be validated. Furthermore, all loads are modelled as constant power loads in DIgSILENT PowerFactory and therefore the same model is assumed in the OPF environment.
As exemplary grid, an existing low voltage feeder in a suburban area in southwestern Germany is modelled. In total, there are 37 points in the grid connecting 14 households.
The simulation is performed in a single time step and only the consumption of the connected households are considered. A power consumption of 6 kVA per household and a power factor of 0.9 is assumed. This comparatively high power demand is assumed to produce high voltage drops in the grid. To cause voltage unbalances, an additional 6 kW load is connected to Phase 1 of two households at the end of the feeder.  Figure 3 and 4, the voltages at several points in the grid can be seen for the OPF presented in this paper and PowerFactory (PF). Figure 3 shows the voltages of Phase 1 and Phase 2. The results for Phase 3 are not shown as they are equal to the results of Phase 2. In Figure 4, the voltage on the neutral line is plotted. In general, the results between the different algorithms are very similar. The maximum deviation is less than 0.0001 p.u. Fig. 4. Voltage on the neutral line IX. EXEMPLARY RESULTS FOR THE UNBALANCED OPTIMAL POWER FLOW ALGORITHM USING A BATTERY STORAGE For the simulations in the following chapter, the same grid data as in Chapter VIII is used. In contrast to Chapter VIII, a simulation for a whole day is performed using a temporal resolution of 15 min. For the household, realistic profiles for a yearly energy consumption of 3000 kWh are used that have been generated using [15]. It is assumed, that the household consumption is balanced. Additionally, each household has charging infrastructure for one-phase charging with 3.7 kW. This admittedly unrealistic assumption is chosen to test unbalanced conditions. The charging profiles for the households are generated using [16]. In Figure 5, the total energy consumption for the test day is shown. The resulting voltages at the end of the feeder are shown in Figure 6. The optimized equipment is a battery storage system with an energy capacity of 100 kWh and power limits of 40 kW is connected at the end of the feeder. The storage system is connected to all three phases and the power flow for each phase can be controlled independently.

A. Ensuring Compliance with the Voltage Limits Using a Battery Storage
In the first test scenario, the maximum voltage drop is restricted to 2 % of the nominal voltage. In Figure 7, the resulting voltages are shown at the end of the feeder. It has to be mentioned that the lowest voltage is not automatically at the end of the feeder. Through discharging (see Figure 8) and charging (see Figure 9) the battery, the voltages in the grid are within the given limits at all time steps.

B. Increasing Self-consumption Using a Battery Storage
In the second test scenario, a solar power system with 5 kWp nominal power is installed at the end of the feeder at Phase 2 in addition to the first scenario. For the profile, measurement data (see Figure 10) from an existing solar power plant in southern Germany is used. It is assumed, that the cost of energy consumption from the utility is 28 ct/kWh, while feed-in tariffs earn only 10 ct/kWh. This is implemented using two instead of one generator per phase at the feeder start. The generator limits were set such that one generator is only importing energy (for 28 ct/kWh) while the second generator is only able to export energy (for 10 ct/kWh). The cost functions are adapted accordingly. Hence, it is economically beneficial to increase the self-consumption in the feeder. As the generation on Phase 2 exceeds the consumption in several time steps, the battery storage is used to store the energy and shift the energy partly to other phases. The resulting energy import is shown in Figure 11. In comparison to Figure 5, the import decreased on Phase 2 significantly and Phase 3 slightly. For the charging Fig. 11. Total power import per phase and discharging processes an efficiency of 90 % is assumed. The battery system is therefore only charged when the solar power generation exceeds the energy consumption on Phase 2 (see Figure 12). The discharging power of the battery system is shown in Figure 13.  In this paper, an optimal power flow algorithm for unbalanced three-phase distribution grids to simulate energy storage systems is presented. The algorithm is validated using a commercial load flow software and two exemplary simulations have been presented.
In the future, the presented algorithm will be expanded to take into account also the negative and zero sequence components of the grid. Additionally, voltage unbalance will be taken into account in the optimization algorithm as an extra grid limit.