Real-Time Active-Reactive Optimal Power Flow with Flexible Operation of Battery Storage Systems

: In this paper, a multi-phase multi-time-scale real-time dynamic active-reactive optimal power ﬂow (RT-DAR-OPF) framework is developed to optimally deal with spontaneous changes in wind power in distribution networks (DNs) with battery storage systems (BSSs). The most challenging issue hereby is that a large-scale ‘dynamic’ (i


Introduction
There is a strong demand for increasing the penetration level of wind energy into distribution networks (DNs).However, a considerable amount of this generation may need to be curtailed due to technical constraints in the network.Battery storage systems (BSSs) can be optimally used to store the energy, decrease the curtailment, and consequently increase economic benefits.However, BSSs introduce dynamic terms to the problem of optimal power flow (OPF).In addition, considering both the active and reactive power capability of the BSSs with flexible operation strategies, as well as maximizing the lifetime of the batteries further increases the complexity of the problem.Furthermore, wind power is intermittent, and therefore the network operator has to quickly update the operation strategies correspondingly.This task should be carried out by a 'real-time' optimization aiming at determining a huge number of mixed-integer decision variables.Therefore, a real-time dynamic active-reactive OPF (RT-DAR-OPF) problem should be solved to ensure not only optimality but also the feasibility of the system operation.
RT-OPF was first introduced in [1] and then attracted the interest of many studies: [2] introduced a real-time strategy based on a linear model predictive control (MPC) for OPF in the presence of BSSs and wind turbines; [3] utilized neural networks for RT-OPF in radial DNs considering uncertain demand and renewable energy generation (REGs); [4] developed two real-time algorithms for constraint satisfaction for OPF; [5] proposed an RT-OPF control approach based on a feedback mechanism; [6] proposed a risk-based RT-OPF considering uncertain REGs; [7] developed a data-driven hourly real-time power dispatch framework; [8] proposed an online gradient algorithm for OPF on radial networks ensuring low computation time; [9] proposed an RT-OPF algorithm based on quasi-Newton methods; [10] proposed an RT-OPF approach considering the minute-to-minute variability of REGs and demand; [11] developed a feedback controller for photovoltaic inverters seeking OPF solutions; [12] extended the research in [11] by improving the convergence properties of the feedback controllers for the case of time-varying ambient and network conditions; [13] developed a prediction-realization RT-OPF framework to deal with intermittent wind power in DNs; [14] extended the framework in [13] by incorporating the reactive power dispatch of wind farms (WFs) as well as overcoming convergence issues [15] by introducing a reconciliation algorithm; [16] utilized graphical processing units to accelerate the OPF computation; [17] presented a real-time optimization method for alleviating contingencies (e.g., line overloads and voltage limit violations) in transmission networks; [18] proposed a multi-stage stochastic optimization for real-time economic dispatch of storage (specifically pumped hydro) resources; [19] proposed a feedback-based RT-OPF algorithm to satisfy technical constraints in real time; and [20] proposed an MPC-based distributed RT-OPF method for interconnected microgrids.For a detailed survey of RT-OPF methods, we refer to [21].It is noted that all the above studies on 'real-time' OPF [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] did not consider the optimal operations of BSSs when optimizing mixed-integer decision variables of the network.
BSSs could play a significant role in the efficient operation of energy networks.They can lead to a decrease in power losses and voltage deviations [22], supplying peak demand [23,24] by storing the energy instead of energy curtailment [25].However, the optimal integration of the BSSs into energy networks still needs further investigations.From an operation point of view, the BSSs can be scheduled in a fixed or flexible manner.For instance, [26] utilized both the active and reactive power capability of BSSs [27,28] in the OPF problem based on a fixed length of charge and discharge periods in a prediction horizon.The reactive power provision of distributed energy resources leads to significant economic (e.g., exporting reactive power to an upstream network [14]) and technical benefits (e.g., management of line losses and voltage regulation [29,30]).In [31], it was proposed that BSSs are operated with a flexible length of charge and discharge periods, but identical operation strategies were obtained for different BSSs so as to reduce the number of mixed-integer decision variables.However, realizing the identical solutions for all BSSs (at different locations) cannot be efficient.Therefore, the work in [31] was extended in [32,33] to find optimal operation strategies for each individual BSS.Recently, [34] proposed a multi-period [35][36][37][38] framework to solve the dynamic OPF problem for distribution networks with BSSs under uncertain renewable energy generation.It is noted that the OPF problems in [26,[31][32][33][35][36][37][38] were not solved in real time and they did not optimize the depth of discharge (DoD) of BSSs while determining flexible operation strategies.Evaluating expended life costs of batteries based on cycle counting only is not conclusive, and therefore, the DoD of BSSs should also be taken into account [39].Flexible operation of BSSs and considering both DoD and the number of charge-discharge cycles in the calculation of expended life costs of batteries lead to more efficient operation of BSSs but highly complex dynamic MINLP AR-OPF, in particular when all mixed-integer decision variables are simultaneously optimized.The incorporation of expended life costs of batteries (in terms of both DoD and number of cycles) into 'real-time' OPF has not been considered yet.Therefore, this study aimed to develop a new approach with the following contributions over the above literature: Energies 2020, 13, 1697 3 of 17 (1) A multi-time-scale dynamic (i.e., with difference equations rather than only algebraic equations) AR-OPF framework is developed to optimally react to the spontaneous changes in wind power and ensure the feasibility of operations in real time when BSSs exist in DNs.(2) The framework offers the possibility of simultaneous optimization of all of the following mixed-integer variables in a prediction horizon: Active-reactive reverse power flow to an upstream network (continuous).
(3) Fully flexible optimal operation strategies for BSSs are determined for the dynamic AR-OPF while minimizing the expended life costs of the BSSs as a function of DoD and the number of charge-discharge cycles.
The remainder of the paper is organized as follows.The formulation of a general stochastic dynamic MINLP optimization problem is provided in Section 2. Section 3 describes the proposed RT-DAR-OPF problem.In Section 4, different modes of operations for BSSs are defined and the dynamic MINLP AR-OPF problem is formulated in detail.The results of a case study and conclusions are provided in Sections 5 and 6, respectively.

Problem Formulation
The aim of the RT-DAR-OPF is to compute optimal operation strategies to be realized to DNs with BSSs under uncertain penetration of wind power.A general formulation of the optimization problem can be expressed as: min u(t),l(t),y(t) f (x(t), u(t), l(t), y(t), ξ(t)) .s.t.
. x(t) = g(x(t), u(t), l(t), y(t), ξ(t)), where f (.) is the objective function to be minimized, x is the vector of state variables, u is the vector of continuous decision variables, l is the vector of discrete decision variables, y is the vector of binary decision variables, ξ is the vector of random variables, and t is time.The objective function is subject to equality and inequality constraints.Here, g(.) denotes dynamic nonlinear model equations, x 0 denotes the initial states at t 0 , x min/max are the lower/upper limits on state variables, u min/max are the lower/upper limits on continuous decision variables, and t f is the final time.Equation ( 1) is a large-scale complex stochastic dynamic MINLP optimization problem, which is difficult to solve.

Real-Time Dynamic AR-OPF Framework
To solve the optimization problem Equation ( 1), we propose a multi-phase multi-time-scale RT-DAR-OPF framework as illustrated in Figure 1.The framework consists of three phases (phases 2 and 3 are adapted from [13,14]) with three different time scales (e.g., T P1 = 24 h, T P2 = 2 min, and T S = 20 s) and 15 steps as follows: (1) Provide hourly forecasted wind power, demand, and price profiles in advance of each prediction horizon T P1 .(2) Solve the corresponding dynamic MINLP AR-OPF problem.In this step, optimal flexible operation strategies for BSSs are computed for the upcoming T P1 (e.g., 24 h, with hourly discretization).
The detailed problem formulation is described in Section 4.
(3) The variables of BSSs computed in phase 1 will be used as fixed input parameters for the second phase.Note that other decision variables will be recomputed in phase 2. ( 4) Provide forecasted values of wind power, demand, and price ahead of each prediction horizon T P2 (e.g., 2 min).Note that the length of the prediction horizon T P2 should depend on the availability of the forecasted data as well as the computation time in step ( 7). ( 5) To describe uncertain wind power, generate N s wind power scenarios for each WF using a continuous bounded stochastic distribution with an identical probability between two adjacent scenarios.For this purpose, N s − 1 intervals are defined for the wind power P w (n w , n s ), n s = 1, . . ., N s , such that: where n w and n s are the indices for WFs and wind power scenarios, respectively.Pr is the probability operator and the scenarios are margins of the defined intervals.To this end, N s wind power scenarios are generated for each WF.Then, N c wind power scenario combinations are formed for each prediction horizon T P2 .The total number of scenario combinations will be: where N w is the total number of WFs.(6) Send the generated N c wind power scenario combinations (obtained in step 5) to the MINLP AR-OPF.(7) Solve the MINLP AR-OPF problems corresponding to each scenario combination for the upcoming T P2 .Note that the optimization problems at this step are not dynamic as the optimal operation strategies of BSSs are already given as input parameters.Since reactive power flow has influence on nodal voltages [40,41], the reactive power dispatch of the WFs can lead to voltage violations, in particular when the wind power fluctuates.For this reason, we use a back-off strategy [14] to satisfy voltage constraints in the RT-OPF.Since the optimization problems in this step are independent, they are solved using parallel computation in order to ensure that the solutions for all the scenario combinations are available within the prediction horizon T P2 .(8) Send the solutions of the MINLP AR-OPF problems (obtained in step 7) as a lookup table to a reconciliation algorithm.(9) Using the reconciliation algorithm, reconcile the lookup table by substituting the un-converged problems with solutions by which the safety of the operations is ensured while minimizing the degree of conservatism.(10) Send the reconciled lookup table for the T P2 to a selection algorithm.(11) Provide the values of wind power measured at each sampling interval T S (e.g., 20 s) to the selection and power factor modification algorithms.
Energies 2020, 13, 1697 5 of 17 (12) The selection algorithm selects a solution strategy based on the measured values of wind power for each sampling interval T S .The selected scenario ensures the safety of the operation with the minimum of the objective function.(13) Send the selected scenario to the power factor modification algorithm.( 14) Modify the power factor of WFs before realizing the solution using the power factor modification algorithm.Due to the possible difference between the measured wind power and the selected scenario, realizing the reactive power dispatch can lead to violations of power factor limits.Therefore, the power factor modification algorithm ensures satisfaction of the power factor constraints.( 15) Send the decision variables to the network at each sampling interval T S .
Energies 2020, 13, x FOR PEER REVIEW 5 of 17 14) Modify the power factor of WFs before realizing the solution using the power factor modification algorithm.Due to the possible difference between the measured wind power and the selected scenario, realizing the reactive power dispatch can lead to violations of power factor limits.Therefore, the power factor modification algorithm ensures satisfaction of the power factor constraints.15) Send the decision variables to the network at each sampling interval S T .
The above 15 steps are repeated for the next 1 P T so that the proposed real-time framework aims to autonomously update the operation strategies according to spontaneous changes of the wind power while optimally managing the operation of BSSs.

Operation Modes of BSSs
It was shown in [39] that the lifetime of a battery is strongly influenced by the DoD and the number of charge-discharge cycles in the prediction horizon.Therefore, in many studies [26,[31][32][33], was limited in order to increase the lifetime.Considering this limitation, three types of operation modes can be defined for BSSs in active-reactive OPF as shown in Figure 2. The above 15 steps are repeated for the next T P1 so that the proposed real-time framework aims to autonomously update the operation strategies according to spontaneous changes of the wind power while optimally managing the operation of BSSs.

Operation Modes of BSSs
It was shown in [39] that the lifetime of a battery is strongly influenced by the DoD and the number of charge-discharge cycles in the prediction horizon.Therefore, in many studies [26,[31][32][33], N cyc was limited in order to increase the lifetime.Considering this limitation, three types of operation modes can be defined for BSSs in active-reactive OPF as shown in Figure 2.  It can be seen that mode 3 allows higher flexibility compared to the other modes of operation.However, it leads to higher complexity as the number of decision variables in the optimization problem increases.Therefore, in this paper, the dynamic AR-OPF is formulated and a solution approach is presented for BSSs to operate in mode 3.For comparison purposes, we also tested our new framework with the other operation modes in the case study.

Detailed Problem Formulation
The dynamic optimization problem in phase 1 of the RT-DAR-OPF framework is formulated as follows: The objective function in Equation ( 4) minimizes the total costs of active and reactive energy 1 F and 2 F imported from an upstream network, and meanwhile minimizes the total expended life costs of the BSSs 3 F .The significance of our proposed dynamic AR-OPF is that all continuous, discrete, Mode 1 [26]: Fixed charge and discharge periods, and a fixed number of charge-discharge cycles in the prediction horizon T P1 .In this operation mode, the decision variables of the BSSs are P ch (i, t), P dis (i, t), and Q b (i, t).
Mode 2 [31][32][33]: Flexible charge and discharge periods, and a fixed number of charge-discharge cycles in the prediction horizon T P1 .In this operation mode, the decision variables of the BSSs are P ch (i, t), P dis (i, t), Q b (i, t), T ch (i, n cyc ), and T dis (i, n cyc ).
Mode 3 (proposed): Flexible charge and discharge periods and optimal number of charge-discharge cycles in the prediction horizon T P1 .In this operation mode, the decision variables of the BSSs are P ch (i, t), P dis (i, t), Q b (i, t), T ch (i, n cyc ), T dis (i, n cyc ), T cyc (i, n cyc ), and N cyc (i).
It can be seen that mode 3 allows higher flexibility compared to the other modes of operation.However, it leads to higher complexity as the number of decision variables in the optimization problem increases.Therefore, in this paper, the dynamic AR-OPF is formulated and a solution approach is presented for BSSs to operate in mode 3.For comparison purposes, we also tested our new framework with the other operation modes in the case study.

Detailed Problem Formulation
The dynamic optimization problem in phase 1 of the RT-DAR-OPF framework is formulated as follows: min The objective function in Equation ( 4) minimizes the total costs of active and reactive energy F 1 and F 2 imported from an upstream network, and meanwhile minimizes the total expended life costs of the BSSs F 3 .The significance of our proposed dynamic AR-OPF is that all continuous, discrete, and binary decision variables are simultaneously optimized for the prediction horizon T P1 .Here, the vector of continuous decision variables u(t) includes curtailment factors of WFs β w (i, t), reactive power dispatch of WFs Q w (i, t), active power charge of BSSs P ch (i, t), active power discharge of BSSs P dis (i, t), and reactive power dispatch of BSSs Q b (i, t).The vector of discrete decision variables I(t) includes slack bus voltage V S (t), charge periods of the BSSs T ch (i, n cyc ), discharge periods of the BSSs T dis (i, n cyc ), length of charge-discharge cycles of batteries T cyc (i, n cyc ), and number of charge-discharge cycles of the BSSs in the prediction horizon N cyc (i).The vector of binary decision variables y(t) includes α(i, t), which represents the status of charge/discharge for each BSS.
Equation ( 4) is subject to the following equality and inequality constraints: where Equations ( 5) and ( 6) are the active and reactive power flow equations at the buses, respectively.
The feeder limits are: S(i, j, t) ≤ S l.max (i, j), i, j ∈ S b ; i j.
The constraints of the curtailment factors of WFs are: and the constraints of the power factors of WFs are:

Equations of BSSs
In this work, with the aid of power conditioning systems (PCSs), the BSSs can provide and absorb both active and reactive power.The active power charge and discharge are constrained to the capacity of the PCSs: Energies 2020, 13, 1697 8 of 17 The reactive power dispatch of the BSSs is constrained to: The apparent power of the BSSs is constrained to: Here, we define a binary decision variable α to avoid charge and discharge at the same time: where α = 1 indicates the charging of the battery while α = 0 denotes the discharging operation.The energy level of a BSS is calculated as follows: where Equation (26) shows the dynamic behavior of a battery with initial states in Equation ( 27).
The ELC of each BSS [39] is a function of the number of cycles in the prediction horizon as well as the average value of the depth of discharge: where: In Equations ( 29) and (33), N cyc (i) can be calculated as follows: where: .

Case Study
In this paper, a 41-bus medium voltage DN [14,26,31,42] is used as a case study to demonstrate the effectiveness of the proposed approach.Two WFs (each with rated power of 10 MW) and two BSSs are located at buses 2 and 16, respectively.The input data for the case study is adapted from [14,39] and given in Table 1 as well as subplots (a)-(c) in Figures 3 and 4. The dynamic MINLP optimization problem in phase 1 is solved using the SBB solver and the optimization problems in phase 2 are solved using the BONMIN solver in GAMS.
Table 1.Data for the case study taken and adapted from [14,39].
Subplots (d)-(l) in Figure 3 show the output of the dynamic MINLP AR-OPF solved in phase 1.The length of the prediction horizon is 24 h with hourly discretization and the maximum number of cycles in the prediction horizon is 4 [39].Based on the forecasted profiles, all the mixed-integer variables are simultaneously solved for the upcoming day.
In Figure 3, subplot (d) shows the optimal values of the slack bus voltage, taking discrete values with the steps of 0.01 pu.Subplot (e) of Figure 3 confirms that using BSSs in the network can lead to a significant reduction of the wind power curtailment.In the same figure, subplots (g)-(i) shows that using the proposed method, the BSSs could operate with dissimilar strategies in the prediction horizon.Subplots (f) and (j) of Figure 3 show that utilizing the reactive power capability of the WFs and BSSs can lead to a huge amount of reactive power generation in the network.Beside the demand and wind power profiles, energy prices play a significant role in determining the optimal operations of BSSs.It means the batteries tend to be charged when the active energy price is low and discharged when it is high.In addition, the BSSs also dispatch reactive power to cover the reactive power demand in the network as well as exporting the surplus amount to the upstream HV network (see Figure 3, subplot (l)).It is noted that in phase 1, only BSSs variables (hourly discretized) are transferred to phase 2 as input parameters.It means the other decision variables are recomputed in phase 2 (with 2-min discretization) and modified in phase 3 (with 20-s discretization) before realization.The results of the realization phase are shown in subplots (d)-(k) in Figure 4.For comparison purposes, we ran the proposed RT-DAR-OPF for three different modes of operation defined in Section 4, with the results shown in Table 2.In modes 1 and 2, the number of cycles per day is fixed to 4, while in the flexible approach, the number of cycles for each BSS is a free variable to be optimized by the solver.
Due to Equations ( 4) and ( 29)-( 36), the optimizer tends to decrease the number of cycles and DoD in order to minimize the expended life costs of the BSSs.However, the effect of the number of cycles on ELC is more significant in our case study as seen in Table 2. Therefore, the total expended life costs of the BSSs are decreased significantly in mode 3 compared to the other operation modes.Moreover, the costs of active and reactive energy at the slack bus are also decreased slightly in mode 3. Altogether, this leads to a a huge reduction in the total costs obtained by using mode 3.In Figure 3, subplot (d) shows the optimal values of the slack bus voltage, taking discrete values with the steps of 0.01 pu.Subplot (e) of Figure 3 confirms that using BSSs in the network can lead to a significant reduction of the wind power curtailment.In the same figure, subplots (g)-(i) shows that using the proposed method, the BSSs could operate with dissimilar strategies in the prediction horizon.Subplots (f) and (j) of Figure 3 show that utilizing the reactive power capability of the WFs and BSSs can lead to a huge amount of reactive power generation in the network.Beside the demand and wind power profiles, energy prices play a significant role in determining the optimal operations

Conclusions
In this paper, a novel multi-time-scale real-time dynamic active-reactive optimal power flow (RT-DAR-OPF) framework was introduced to deal with fast-changing wind power in the presence of battery storage systems (BSSs).The framework consists of three phases: In the first phase, a dynamic mixed-integer nonlinear programming problem is solved to simultaneously determine the optimal operation strategies of the BSSs for the upcoming day.In the second phase, wind power scenarios are generated based on the forecasted wind power values for short prediction horizons (e.g., 2 min) and then the AR-OPF problems corresponding to the scenarios are solved in parallel.The results are saved as a lookup table from which one solution is selected based on the actual values of wind power in a very short sampling time (e.g., 20 s) in the third phase.The solution is then modified to ensure satisfaction of the constraints.The operation strategies obtained by the proposed fully flexible optimal operation strategies of BSSs show significant advantages over the results of the methods with fixed operation strategies of BSSs.This is mostly due to the reduction of the expended life costs of the batteries in the proposed method.Since the framework safeguards the feasibility and optimality of the operations in real time, it could be used for real applications in the future.

Figure 3 .
Figure 3. (a) Total active and reactive power demand; (b) Energy prices; (c) Wind power; (d) Slack bus voltage; (e) Curtailment factors of WFs; (f) Capacitive reactive power of WFs; (g) Binary variables for charge/discharge of BSSs; (h),(i) Active power charge and discharge of BSSs, respectively; (j) Capacitive reactive power of BSSs; (k) Energy levels in BSSs; (l) Active and reactive power at the slack bus.

Figure 3 .
Figure 3. (a) Total active and reactive power demand; (b) Energy prices; (c) Wind power; (d) Slack bus voltage; (e) Curtailment factors of WFs; (f) Capacitive reactive power of WFs; (g) Binary variables for charge/discharge of BSSs; (h),(i) Active power charge and discharge of BSSs, respectively; (j) Capacitive reactive power of BSSs; (k) Energy levels in BSSs; (l) Active and reactive power at the slack bus.

Figure 4 .
Figure 4. (a) Total active and reactive power demand; (b),(c) Actual wind power of the first and second WF; (d),(e) Curtailment factors of the first and second WF, respectively; (f) Capacitive reactive power dispatch of the first WF; (g) Capacitive reactive power dispatch of the second WF; (h) Slack bus voltage; (i) Active power at the slack bus; (j) Reactive power at the slack bus; (k) Total costs of active and reactive energy at the slack bus.

Figure 4 .
Figure 4. (a) Total active and reactive power demand; (b),(c) Actual wind power of the first and second WF; (d),(e) Curtailment factors of the first and second WF, respectively; (f) Capacitive reactive power dispatch of the first WF; (g) Capacitive reactive power dispatch of the second WF; (h) Slack bus voltage; (i) Active power at the slack bus; (j) Reactive power at the slack bus; (k) Total costs of active and reactive energy at the slack bus.

Table 2 .
Comparison of the RT-DAR-OPF with the proposed flexible operation strategy to the results obtained by mode 1 and mode 2 for one day.Operation Mode N cyc (1) N cyc (2)