Incentive Price-Based Demand Response in Active Distribution Grids

: Integration of PV power generation systems at distribution grids, especially at low-voltage (LV) grids, brings in operational challenges for distribution system operators (DSOs). These challenges include grid over-voltages and overloading of cables during peak PV power production. Battery energy storage systems (BESS) are being installed alongside PV systems by customers for smart home energy management. This paper investigates the utilization of those BESS by DSOs for maintaining the grid voltages within limits. In this context, an incentive price based demand response (IDR) method is proposed for indirect control of charging/discharging power of the BESS according to the grid voltage conditions. It is shown that the proposed IDR method, which relies on a distributed computing application, is able to maintain the grid voltages within limits. The advantage of the proposed distributed implementation is that the DSOs can compute and communicate the incentive prices thereby encouraging customers to actively participate in the demand response program. An iterative distributed algorithm is used to compute the incentive prices of individual BESS to minimize the costs of net power consumption of the customer. The proposed IDR method is tested by conducting simulation studies on the model of a Danish LV grid for few study cases. The simulation results show that by using the proposed method for the control of BESS, node voltages are maintained within limits as well as the costs of net power consumption of BESS owners are minimized. achieved by employing an iterative distributed algorithm that computes the incentive prices iteratively. Simulation results show that the costs of electricity power consumption for the customers with BESS is less if the proposed DR method is used. The studies in this paper have not considered the other objectives of the LVGC such as power loss minimization. However, they can be included in the proposed method after the inclusion of the corresponding objective functions and constraints. Future work includes extending the proposed incentive price based DR for EVs, electric boilers and refrigeration systems.


Introduction
At present, the medium-voltage (MV) and low-voltage (LV) power distribution grids are being integrated with large amounts of renewable energy sources (RES) such as solar PV and wind power plants to accelerate the green energy transition [1]. New types of active loads such as electric vehicles and heat pumps and battery energy storage systems (BESS) collectively called distributed energy resources (DERs) are being increasingly used by customers at LV distribution grids [2]. Because of the above reasons, operation of the distribution grids in the near future will be challenging for the distribution system operators (DSOs). Some of the operational challenges are insufficient hosting capacity of the grid for integrating more RES, overloading of substation transformers and cables, under-voltage and over-voltage events [3][4][5]. In LV grids, the penetration of kilo Watts level PV systems is increasing and is expected to be high enough to cause reverse power flows during peak PV power production. As residential loads are generally predominant in LV grids, the power demand will have peaks in the morning and evening and it may be less in the noon when the PV generation is at its peak. This could lead to over-voltage events during midday and under-voltage events in the morning and early evening. It is to be noted that as of today, PV inverters are operated at unity pf to maximize the PV power production. The grid codes mandate reactive power capability of PV inverters; e.g., in [6] it is stated that PVs with ratings above 11 kW should support voltage regulation in distribution grids. However, the provision of reactive power reduces the active power producing capacity, controller (LVGC) in this paper which is operated by the DSO to control all DERs in their LV grid [13]. In hierarchical distributed control method, the LVGC acts as a coordinator and communicates with the DER controllers to compute the optimum setpoints of DERs. A centralized DR control is typically applicable for grids under normal conditions in which DSOs want to utilize DERs to increase the grid energy efficiency. The centralized control method is not scalable, hence it is not feasible for large distribution grids with thousands of DERs. The centralized method is also not reliable due to its dependability on a single central controller. The above demerits are overcome in the distributed control method as this method involves cooperation among controllers and the failure of few controllers will not cause a large disruption. The hierarchical distributed control can be used for addressing the DSO's objectives such as to mitigate the over-voltages and for meeting the customers' requirement to minimize their electricity bills [14,15].
The focus of this paper is to apply a novel DR method for two objectives, which are voltage regulation of LV grids and maximization of the self-power consumption of customers (thereby minimizing their electricity bills) by optimum utilization of their battery storage systems. Distributed model predictive control (DMPC) method can be employed for realizing a distributed DR scheme as discussed in [15] to make use of the forecasted power demand for calculating the optimum control setpoints to DERs. This is the motivation for using a DMPC method for the distributed implementation of the DR scheme in this paper. In order to dynamically compute the incentive prices that can solve the voltage regulation issues price-based distributed optimization methods [16][17][18] are to be considered. An incentive-based distributed optimization algorithm is proposed in [18] for voltage regulation by optimum control of the PV power production in LV grids. This algorithm is proposed to use primal-dual projected gradient method to compute the incentive price signals of PV systems. However, the disadvantages of this approach are (i) the PV power will be curtailed to bring the grid voltage within limits, (ii) forecasted load and generation information are not utilized for computing the incentive price signals. The above demerits can be overcome as follows (i) charging of the BESS when there is an over-voltage instead of curtailing the PV power, and (ii) use of forecasted power demand to solve the optimization problem over a prediction horizon, thereby the optimizing the utilization of BESS over a future time period. The benefits of the proposed approach will be increase in the hosting capacity of the grid, reduced operating costs for the DSO, and reduced electricity bills for the customers. A profit-sharing model or a detailed economic analysis among the players are not considered and they are beyond the scope of this paper. A limitation of the proposed method is the requirement of communication infrastructure among the LVGC and BESS controllers and the need for exchange of signals at each time-step to compute the incentive prices. However, it is assumed that in the near future, with the rollout of smart meters and other smart grid technologies the above communication requirement could be satisfied.
The following are the assumptions of this work.
(A1) Hourly load forecasts and PV generation forecasts are available that can be used by both DSO and BESS controllers. (A2) The LVGC and DER controllers are able to solve the proposed distributed optimization routine at each time-step (10 min) by exchanging information such as dual variables of the optimization problem and the incentive price signals over the communication network.
It is to be noted that short-term forecasts of loads and generation are important to realize a model predictive control algorithm, which is stated in assumption A1. As per assumption A2, it is assumed that the DSO install and operate the LVGC, which may be located in the LV distribution substation. Similarly, the DERs owned and operated by the customers are assumed to have the computational power and communication capability with the LVGC. In this paper, the challenge of DSOs to compute incentive price signals to customers owning flexible resources such as BESS is addressed. The major contributions of this paper are provided below.
1. Primal-dual projected gradient-based distributed optimization algorithm [18] is proposed to calculate the incentive price signals by the DSO to regulate the node voltages within limits. 2. Model predictive control algorithm for control of individual BESS to minimize the electricity costs of the corresponding customer taking into account the incentive price signals from the DSO.
It is to be noted that compared to [18], a more accurate method for estimating the node voltages is employed in this paper. Although Ref. [13] also employs a model predictive control algorithm for control of individual BESS, the main difference is the distributed implementation of the above algorithm instead of a centralized control approach.
This paper is organized as follows. A review and comparison of the demand response methods for the control of DERs are made in Section 2. In Section 3, linear models of the LV grid and the BESS are derived and the simulation setup used for conducting simulation studies in this paper are described. Section 4 elaborates the formulation of the proposed incentive price based demand response control method of DERs. Three simulation case studies to numerically evaluate the proposed IDR method are reported in Sections 5 and 6 concludes the paper with remarks and future works.

Demand Response Methods for Control of DERs: A Review and Comparison
In this section, a review of demand response methods for close to real-time of DERs as well as a comparison of their merits and demerits are presented. Demand response of DERs can be implemented either in a centralized or distributed manner [12]. Both these methods are reviewed as follows. Figure 1 shows a demand response method based on solving a centralized optimization problem [13]. This methodology employs a central controller called low-voltage grid controller (LVGC) located either at the DSO control center or at the secondary substation. This controller computes power setpoints to each DER at every time-step (typically 10-15 min) based on the received feedback from the grid, states of the DERs (e.g., the state of charge of the BESS) and the forecasted power values. Each DER has a local controller called DER controller, which receives the power setpoints from the LVGC and does the local control of the DER. The distribution system state estimation (DSSE) algorithm takes fewer measurements from the grid at critical nodes, and uses pseudo-measurements at rest of the nodes to estimate the grid states i.e., the node voltages as well as the powers at all nodes [19,20]. A forecasting algorithm provides hourly forecasted values of the PV generation and the load consumption at all nodes for the next few hours. The typical control objectives of the centralized optimization problem are minimizing the grid power losses, voltage fluctuations over time, maximizing the self-power consumption of PV power generation etc. [21].  As seen in Figure 1, the centralized controller needs to know the grid states (node voltages and powers), and internal states of the DERs (e.g., state-of-charges of the BESS and room temperatures of the home heating systems controlled by heat pumps). This method may be useful for demand response applications in small and medium-scale distribution grids with small number of DERs.

Centralized Demand Response Control of DERs
The demerits of the above centralized control method are as follows.
• Data privacy of the customers is not ensured as the states of DERs contain sensitive information. • The centralized control method will have issues with respect to scalability and reliability and therefore cannot be applied to larger distribution grids with thousands of DERs. • The customers are not empowered to be active and dynamic when they participate in demand response programs. • The communication overhead will be high as this method involves communication between DERs and the DSO control center.

Distributed Demand Response Control Methods of DERs
To overcome the drawbacks in the centralized DR control method, distributed versions of the centralized methods are proposed in the literature [22][23][24][25]. In distributed optimization based DR methods, the centralized optimization problem will be decomposed into two parts assuming the original problem is separable. Two different optimization problems are solved iteratively by the LVGC and DER controllers independently by exchanging of the so-called dual variables and then combining the results to reach a solution closer to that of the centralized problem. The benefits of the distributed optimization approach are the following.

•
Both the DSO and the customers can collectively work for an economical demand response. • Computational requirements for the LVGC will be less and the method will be scalable and extensible for various types of DERs. • Customers will have more freedom in terms of their availability, quality of service etc., during their participation in demand response programs.
There are two variants of the distributed optimization based DR methods which are either based on energy commitment [26,27] or calculation of incentive prices for indirect control of DERs [16][17][18].
1. Energy commitment based DR: The LVGC of the DSO and DER controllers of the customers mutually calculate the optimal amount of flexible energy require for a future time period such that the overall costs of operation of the network as well as electricity bills of the customers are minimum. 2. Incentive prices based DR: The LVGC of the DSO and DER controllers of the customers mutually calculate how much incentives in electricity prices are to be provided by the DSO such that the overall costs of operation of the network as well as electricity bills of the customers are minimized.
It is to be noted that both these methods require a certain number of iterations between the LVGC and the DER controllers to find the optimal flexibility. The difference between these two methods is that the DSO computes directly the flexible energy in the first approach while the incentive price signals are calculated in the second approach.
Though the first approach based on energy commitment mentioned above is better than the centralized optimization method, it has the following drawbacks.

•
DSOs could not motivate the customers to provide more flexibility than they are prepared to offer as the iterations involve only around requested flexible power from DSO and available flexibility from DERs. • Customers may be motivated to provide more flexibility, if the DSO uses incentive prices as the variables to negotiate with the customers.
This paper proposes an incentive price-based demand response method based on the second approach to overcome the above disadvantages. The schematic diagram of the proposed IDR method is shown in Figure 2 and an explanation of its execution is provided below.

•
An iterative distributed algorithm computes the incentive price signals for each BESS at each time-step. The objective of this algorithm is to minimize the costs of net power consumption over time in the nodes where BESS are connected and to maintain the grid voltages within allowed limits. These iterations take into the consideration of the BESS constraints such as SOC limits, power limits and the grid constraints such as voltage limits while calculating the incentive price signals. • Once the iteration converges within a specified tolerance, the computed incentive prices are communicated to the BESS. The BESS controllers use the incentive signals to decide the charging/discharging power of the BESS until the next time-step. • The incentive price signals will be zero when there are no grid voltage issues, hence the algorithm is iterated only once. In this case, no flexibility is requested from the BESS controller. When incentive price signals are nonzero, the BESS controllers are requested to provide flexibility based on their available power rating and the state of charge.  The proposed IDR method has the merits of directly computing the incentive prices and it is superior to the centralized method because customers are incentivized to participate in the demand response program for their energy cost savings. An iterative distributed algorithm based on [18] is proposed to calculate the incentive price signals that can be communicated to the DERs for demand response.

Mathematical Modeling and Simulation Setup
In this section, modeling of a balanced three-phase LV distribution grid of radial topology is presented. Let N + 1 be the total number of nodes of the LV grid including a slack node which represents the secondary side of the LV distribution transformer. Let the slack node be denoted as 0 and the nodes at which the loads are connected be denoted by the set N := {1, . . . , N}. As the LV grid is considered to be of radial topology, there will be N number of branches. Let the branches of this grid be represented by the set L with N × N elements which are arbitrarily numbered and are oriented outwards from the slack node.
It is to be noted that the line resistances in distribution grids are higher than the reactances and needs to be included in the model [28]. The impedance of each branch (i, j) ∈ L is denoted as z ij = r ij + jx ij , where r ij and x ij are the branch resistance and reactance respectively in pu. The shunt admittances at each node are not included in the model as their values are negligible due to short line lengths. Let R db = diag(r ij ) (i,j)∈L and X db = diag(x ij ) (i,j)∈L be the diagonal matrices with branch resistance and reactance values respectively.
Let us consider the slack node voltage to be the reference, i.e., arg(v 0 ) = 0. Let us denote the node-to-ground voltage and current injected at the N nodes by the complex phasorsv n ,ī n and be collected in the setv := [v 1 , · · · ,v N ] T andī := ī 1 , · · · ,ī N T respectively. Let the vectors p and q contain the values of net active power and reactive power demand (power consumption minus generation) at each node and are of size N.

Linear Model of the LV Grid
A linear model of the LV grid based on [13] is derived below to express the node voltages as a linear function of node active and reactive powers. The branch currents i b := [ī b1 , · · · ,ī bN ] can be computed from the bus currentsī as below where the binary matrix M r is the reduced version of the bus injection to branch current matrix with the column pertaining to the slack bus current removed [29]. As per Ohm's law, the branch voltage drops can be computed from the bus current injections as Then each node voltage phasor can be written as the difference between the slack node voltage v 0 and the voltage drop on the line between the corresponding node and the slack node. In vector format, it is written as follows.
where the vector v s := [v 0 , · · · , v 0 ] T of size N. From (1), the following expression can be written.
Let us assume that voltage magnitudes of the N buses are either measured or estimated at each sample time, denoted as v m = [v m1 , · · · , v mN ] T and collected in the diagonal matrix V dm = diag(v m ). Then the approximate expression for branch currents from (1) can be written in terms of the bus active and reactive powers as below.
where the bus current injections are calculated using the approximationī = V −1 dm (p + jq) * . The real part of the node voltages from (3) can be expressed as below.
Using (4), the above expression can be written in terms of bus active and reactive powers as below.
In LV grids, the voltage phase-angles are small, hence the imaginary part of the node voltages from (3) can be neglected [28]. The errors introduced by this approximation are very much negligible (the maximum error will be less than 0.4%). Based on this assumption, (6) can be written as Let us define the constants The above expression for bus voltage magnitudes is linear with respect to the bus active and reactive powers. The matrices B vp and B vq are updated at each sample time based on the latest voltage measurements/estimates.
In [13], it is shown that the above formulation for estimating the node voltages is more accurate than the state-of-the-art methods such as LinDistFlow [30].

Linear Model of the BESS
In this paper, a simplified discrete-time linear model of the BESS based on [31] is considered. The state-of-charge of the battery denoted by x f is related to the battery discharging power p f as follows.
where κ is the time instant, η f is the battery efficiency, E f is the battery capacity and τ s is the sampling time.

Model Predictive Control Algorithm of the BESS
This paper assumes that each BESS is controlled to maximize the self-consumption of the power at the node to which it is connected. A residential load and a PV is assumed to be connected along with each BESS. A sophisticated controller is assumed to be responsible for charging/discharging control of each BESS. Let N f ⊆ N be the nodes where, the BESS are connected. Let the total number of BESS be denoted as N f , hence N f = {1, · · · , N f }. Let the control horizon be denoted as N c which is assumed to be equal to the prediction horizon N p and let the set N c = {1, · · · , N c } contains the discrete-time instants of the control horizon. The net power consumption at each node and at each time instant κ is calculated from the difference between net power demand which is the sum of power consumption of the load (p d ) and the BESS power (p f ), and the PV power generation (p pv ) as provided in (10).
It is to be noted that a typical load profile for the power consumption and a typical power generation profile for the PV on a clear sky day are considered in this paper.
The MPC of the BESS involves solving the below optimization problem to compute the active power setpoints of each BESS denoted as p f .
subject to (10) and the following constraints, where κ is the time instant. The above objective function tries to compute the power from BESS that will minimize the costs of power consumption at each node over the control horizon, thereby maximizing the self-power consumption over time.
The constraint that specifies the state evolution of the BESS over a control horizon using the simplified linear model expressed in (9) is provided below.
. The limits of the battery SOC are stated in the below constraint.
The values of the BESS parameters are chosen to be E f = 20 kWh, η f = 0.91 pu and τ s = 10 min. The active power limits of the BESS is enforced in the below iteration. As we assume the BESS to be operating in unity pf mode, the maximum limit of the BESS's active power is equal to its apparent power rating.

Proposed IDR Method for Control of the BESS
This section presents the formulation of the proposed IDR, which employs a distributed optimization problem. Unlike in the centralized implementation, the LVGC calculates incentive price signals for individual BESS over a prediction horizon such that the objectives of both LVGC and BESS controllers are met. The optimization problem of LVGC is presented first followed by that of the BESS controllers.

Proposed Optimization Problem of the LVGC
The typical objectives of the LVGC are minimization of the grid active power losses, optimum use of DERs and optimum active and reactive power drawn/exported to the MV grid subjected to the constraints of the LV grid such as node voltage limits, branch current limits, and power limits of the DERs etc. In this paper, the discussion is limited to the optimum flexible power to be provided by DERs for maximizing the self-power consumption at each node by providing incentive price signals to each DER. Although the same formulation can be extended to include other objectives as well, it is out of scope of this paper and will be taken into consideration in our future works.
Let us assume that the DSO can provide an incentive price signal α κ f i to each customer who owns the BESS denoted by i at each time instant κ. Each customer may react to the incentive price signal differently and the DSO cannot influence the customer's response to provide flexibility. As a worst case, the customer may ignore the incentive price and will not provide any flexibility. Alternatively, their response could be a linear variation of BESS power consumption with respect to the incentive price signal over a certain range and beyond that there will not be any increase in power consumption due to limitation in the power rating of the BESS. Let us formulate the above dependency by an unknown function p κ f i = f i (α κ f i ). A simplistic assumption, as mentioned earlier, is that it could be a ramp function with bounds on both ends as stated in [18], but this assumption will make the function non-affine. The above unknown function could also be nonlinear depends on configuration of the DER controller's response to the incentive price signal.
The objective of the LVGC is defined as the minimization of the costs of net power consumption p n at each node n over the control horizon, which indirectly means maximization of the self-consumption of power at each node, by the following expression subject to the following constraints in which κ ∈ N c . The matrix Q κ is a diagonal matrix containing the costs of electricity power consumption at the time instant κ. The constraint specifying the node voltages as a linear function of node active and reactive powers based on the derivation in Section 3.1 is provided below.
The minimum and maximum limits of node voltages are specified in the below constraint whose typical values are 0.9 and 1.1 pu respectively.
The power limits of each DER is forced in the below constraint.
The below constraint specifies the dependency of the power consumption of the ith DER on the incentive price signal α i provided by the DSO by an unknown function f i at each time instant κ.
This function need not be an affine function because the incentive price calculated for each BESS need not have a linear relationship with its power consumption. This property makes the optimization problem of the LVGC formulated in (12) to be non-convex.
A convex relaxation of the above optimization problem is derived in [18] in which the constraint (12e) is removed and the modified formulation is provided below.
subject to the following constraints in which κ ∈ N c . As the reactive power capability of BESS is not utilized in the proposed method, its value is set to zero in (12b).
Instead of the nonconvex constraint (12e), the below constraint is used which simply states the power of the BESS at any time should within the feasible values collected in the set P f and for a BESS it will be its maximum and minimum possible power values.
As per the technique used for the convex relaxation of (12e) in [18], the incentive price signals for the BESS can be computed from the below expression. The detailed mathematical derivation can be found in [18].
where L bs is a binary matrix that maps the location of each BESS to the grid nodes, µ v min and µ v max are the dual variables of the node voltages.

Proposed Optimization Problem of the BESS
This paper proposes the below optimization problem to be solved by each BESS in order to find their charging/discharging powers at the present and future time instants, i.e., over the prediction horizon N p . The optimization problem for control of DERs defined in (11) is modified to include incentive price signals in the objective function as shown in (14).

Distributed Optimization Based on Incentive Prices
The optimization problems of the LVGC and BESS controllers defined in (13) and (14) respectively can be solved by the below iterative distributed algorithm. The proof of the following can be found in [18] and it is omitted here.
Update of power setpoints by each flexible resource are done by the below expression, where k denotes the iteration number.
Update of dual variables of the node voltages are made according to the below equation.
where [ ] + denotes the projection operation onto the set R N + , φ > 0 is a regularization factor, and ε 1 , ε 2 > 0 are predefined parameters to prescribe the step sizes for the dual variables update. The node voltages are updated with the latest value of flexible powers from BESS as follows The incentive price of flexible powers from BESS are updated from the values of updated dual variables.
The steps for the iterative algorithm to calculate incentive price signals is provided below. As shown in Algorithm 1, the LVGC and BESS controllers exchange information (incentive price signals and dual variables) until the iterations converge within a specified tolerance.

Algorithm 1
Distributed iterative algorithm to compute incentive prices [18] Initialize the variables: BESS active powers p f , dual variables µ min and µ max and the incentive prices for BESS active powers α f . for{m = 1,2, · · · (repeat until convergence)} do

Simulation Studies
In this section, the simulation studies done for a numerical evaluation of the proposed IDR method formulated in Section 4 are presented. A simplified model of a LV distribution grid from Askov, Denmark with 17 nodes connected with residential customers, PVs and BESS is used for the case studies done in this paper. The schematic diagram of this LV grid is shown in Figure 3. The cable parameters and the sizes of PV systems and BESS connected at each node of this grid can be found in [13]. As of today, the customers connected to this grid have no BESS systems. The PV systems are operated with unity pf and there are no voltage problems experienced in this grid. More PVs as well as BESS alongside some of the PVs are considered for the simulation studies in this paper based on a future operational scenario of this grid. The penetration level of PVs calculated from the ratio of the rated PV generation to the rated load power is around 95% and the total ratings of all BESS is 62% of the LV substation transformer rating. The battery energy capacity is set to be 20 kWh for all the BESS. The purpose of the above addition is to evaluate the grid for potential over-voltage issues without and with the proposed demand response control of the BESS. A typical daily profile of PV power generation and load consumption with 10 min time resolution is shown in Figure 4. It is to be noted that the load active and reactive power profiles and the solar irradiance profile used in this work are provided by the DSO, Norlys. All the case studies done in this paper are done using the Optimization Toolbox™ of Matlab/Simulink™ on an Intel ® Core-i7 2.59-GHz personal computer.

Case 1: Simulation without BESS Connected to the LV Grid
The aim of this case study is to analyze the voltage profile of the LV grid in Figure 3 with typical PV power generation on a sunny day. A simulation is performed without any BESS connected to the LV grid and it will be a base case used for comparison of the results with that of the case with BESS. A typical daily profile of PV power generation and load consumption shown in Figure 4 is used for this study. An increased load consumption especially during evening is used in this simulation considering a future operational scenario. The PV systems are operated with unity pf, which is typically the case in few European countries such as Denmark.
The plot of voltages at the end nodes (R7, R13 and R17) are shown in Figure 5. During the peak PV production from 10:00 to 15:30 h, the node voltages are higher than 1.1 pu which is the allowed limit as per EN50160 standard [32]. The nodes voltages are below the lower limit of 0.9 pu from 18:30 to 21:00 due to high load consumption. It is to be noted that the voltage violations are noticed in this simulation study because of an increase in the penetration of PV generation to 95% and the load consumption. However, the real LV grid, as of today, does not experience any voltage problems.

Case 2: Simulation of BESS Control without Proposed IDR Method
In this simulation case, BESS are connected to the 10 nodes of the LV grid as specified in Figure 3. The BESS are controlled to maximize the self-power consumption at each node i.e., the excess power generated is stored in the BESS during peak PV power production and the BESS is discharged during off-peak hours. The power setpoints of the BESS are computed by solving the optimization problem defined in (11) whose objective is to minimize the cost of net power consumption at each node. The prediction horizon is set to 4 h (N p = 4). The choice of prediction horizon depends on many factors such as the computational burden and memory storages of the LVGC and BESS controllers, duration of the available forecasted powers, and the communication bandwidth. Although choosing a higher prediction horizon improves the cost savings slightly, in this study a 4 h prediction horizon which is fair enough to address the over-voltage issues due to PV power production. The lower and upper limits of SOCs of the BESS are set to be 0.1 and 0.9 respectively.
The simulation results can be seen in Figure 6. The plot (a) shows the power setpoints of BESS every 10 min. The SOC of the BESS at end nodes (R7, R13 and R17) are shown in the plot (b) while the node voltage profiles are shown in the plot (c). It can be seen that the BESS is fully utilized for minimizing the power export to the grid and their SOCs reach the upper limit of 0.9 pu. However, the node voltages are slightly above the limit of 1.1 pu during noon and it is below the limit of 0.9 pu in the evening.  The BESS charging and discharging routine in this case study are independent of the grid voltages, hence there are over-voltage and under-voltage events similar to case 1, although the magnitude and duration of these voltage violation are less. If the storage capacity of the BESS is increased, the voltage violations can be indirectly managed better. However, the above solution is expensive and there is no guarantee that the node voltages are always maintained within the limits.

Case 3: Simulation of BESS Control with the Proposed IDR Method
In case 2, the BESS are controlled to maximize the self-power consumption, which does not take into consideration of the grid voltages. As a result, the end nodes of the LV grid has over-voltages during noon and a slight under-voltage event in the evening. In this case study, the proposed IDR method which uses the distributed optimization technique as per (15) is used to compute the incentive prices of the BESS. Similar to case 2, the value of N p is set to 4.0. At each time-step, the incentive price for the operation of each BESS is calculated using the procedure provided in Algorithm 1.
The simulation results are shown in Figure 7. The first plot shows the incentive prices for the BESS at the end nodes. From the results of case 1 shown in Figure 5, it can be seen that over-voltages at end nodes occur after 10:00. Due to this fact, the calculated incentive prices for BESS are zero until that time in this case study. It is to be noted that a positive incentive price indicates reduction in electricity price for the customer, thus encouraging the customer to charge the BESS. On the other hand, a negative incentive price denotes increase in electricity price to the customer encouraging them to reduce the net electricity consumption, which could be achieved by discharging the BESS, if possible. Due to the calculated positive incentive prices during peak power production, the voltages are kept within the limits. The under-voltage event, which would have occurred during evening peak load period, is also avoided by generation of negative incentive price signals to the BESS during that period.
The computation time for the execution of case 2 at each time-step is 0.3 s, while for case 3 it is about 0.6 to 3.26 s on a desktop computer representative of current standards. The execution time is higher in case 3 because of the iterations involved in the computation of incentive price signals. It is assumed in the simulation case studies that the uncertainties in the system parameters, measurements etc., are not high to cause a major difference in the results obtained. Also, the proposed algorithm is executed once in 10 min interval taking into account the updated measurements from the LV grid, hence any deviations in the demand response control caused by uncertainties will be corrected in the next time-step. The results of three simulation case studies provided in Table 1 lead to the following conclusions.

•
The proposed method may be more advantageous to the DSOs than the time-of-use tariffs in which the electricity prices are normally predetermined and fixed. However, in the proposed IDR method, the incentive prices are dynamically calculated based on the grid operating conditions. Hence using the proposed method, the DSO can minimize their operating costs and increase the hosting capacity of the grid. • From the customers' point of view, the proposed method is able to make use of the incentives from the DSO to optimally use their DERs and reduce their electricity bills.

Conclusions
This paper proposes a demand response method in which the incentive prices are computed and communicated to each DER close to real-time by which the voltage regulation objective of the DSO and electricity cost minimization of the customer are both fulfilled. The proposed method is verified through simulation studies in which it is applied for controlling the charging/discharging pattern of a group of BESS connected to the simplified model of a real-life LV grid with PV systems. In this paper, the optimization problems of the LVGC and BESS controllers are formulated over a receding horizon using the forecasted powers. The distributed implementation of the above-mentioned optimization problems is achieved by employing an iterative distributed algorithm that computes the incentive prices iteratively. Simulation results show that the costs of electricity power consumption for the customers with BESS is less if the proposed DR method is used. The studies in this paper have not considered the other objectives of the LVGC such as power loss minimization. However, they can be included in the proposed method after the inclusion of the corresponding objective functions and constraints. Future work includes extending the proposed incentive price based DR for EVs, electric boilers and refrigeration systems.
Author Contributions: K.N., the main author of this article, prepared the draft. J.R.P. and B.B.-J. the co-authors, have provided technical inputs, reviewed the draft and supervised the work. B.B.-J. was also responsible for the administration and funding acquisition of the project "DECODE". All authors have read and agreed to the published version of the manuscript.
Funding: This research work was funded by Energinet, the Danish transmission system operator, under PSO-ForskEL programme, grant no. 12414, for the project "DECODE".