Predictive Control of Flexible Resources for Demand Response in Active Distribution Networks

In this paper, a model-based predictive control method is proposed for utilization of flexible resources such as battery energy storage systems and heating systems effectively to provide demand response in low-voltage distribution networks with solar photovoltaic. The contributions of this paper are twofold. First, a linear power flow method based on relaxation of branch power losses applicable to radial distribution networks is proposed and formulated. Second, a flexible resources controller that solves a multi-objective linear optimization problem in receding-horizon fashion is formulated taking into account system states, forecasts of generation, and loads. Using the proposed control algorithm, flexibility from network resources can be utilized for low-voltage network management with assurance of quality of service to the customers. Simulations are conducted for summer and winter cases on a simplified Danish low-voltage network using Matlab/Simulink to study the performance of the proposed control method. Compared to the methods in state of the art, the proposed linear power flow method is proven to be accurate for the calculation of network power flows. Simulation results also show that proposed flexible resources controller can meet the network control objectives while satisfying the network constraints and operation limits of the flexible resources.


I. INTRODUCTION
T HE modern power distribution system is undergoing a phenomenal change due to green energy resources and active loads with flexible power consumption being integrated to the medium-voltage (MV) and low-voltage (LV) active distribution networks (ADN) [1].Utilization of the flexible resources present in the network for demand response (DR) offers many advantages including accommodation of high penetration of distributed energy resources (DER) such as solar photo-voltaic (PV) and wind power generators, postponing or avoiding the network reinforcement and prevention of active power curtailment from DER [2]- [4].Network congestions which may be caused due to high penetration of DER in distribution networks can be managed by employing a suitable control method to monitor and control flexible resources [5].
In this paper, the term flexibility of a resource refers to the time shift and/or magnitude change in power consumption in response to an external control signal.The architecture and control methodology for achieving cost-effective operation of the distribution network by utilizing the available flexibility is an active research topic and is also the focus of this paper.
Hierarchical control structures are proposed in the literature [6], [7] for cost-effective operation of the active distribution networks.In [6], an on-load tap changing transformer at an HV/MV substation and DERs at the LV network are used to provide flexibility.An optimal method to control a group of DERs as virtual power plants is proposed in [8].In [7], electric vehicles are utilized as flexible resources for energy management.However, to accommodate more renewable energy from DER without overvoltages and overloading of lines, it is essential to utilize the flexible resources present in the LV network.Based on the analyses in [3] and [9], in cold countries such as Denmark, there is a good potential to utilize heat pumps (HP) at residential buildings of LV networks to provide short-term flexibility.Also, battery energy storage systems (BESS) installed along with PV systems can be utilized as a flexible resource for grid voltage support, peak shaving and energy management [10].While utilizing the flexibility from above resources for demand response, the control method should assure the quality of service for the customers [11].
In this paper a multi-objective model predictive control (MPC) method is proposed for utilizing the flexibility from BESS and HP in LV distribution networks.The objective is to control BESS and HP based on their linear models within the operational constraints, by solving an open-loop optimal control problem on-line over a finite horizon.The proposed control method make use of the estimated states of the systems and forecasts of the output disturbances to compute the current control action.It is assumed that observability modules are available to provide system states and forecasts information.The main motivation to use MPC method in this work is to make use of this available network information to predict the future power flows of the network and take appropriate control actions at the present time instant.To include the network power flow model in the above formulation, linearization of the ac power flow equations is necessary.Various methods to linearize the ac power flow equations are reported in the literature [8], [12] and [13].
A survey of control methods for optimal power flow in power networks can be found in [14] and the references therein.Economic model predictive control (MPC) formulated using linear programming (LP) technique is applied in [9] to find the optimal operation of HP to reduce electricity costs.However, the influence of network power flows and provision of demand response by HP are not investigated.In [10], an LP method is used to find the charging profile of the BESS to maximize the power absorption from solar PV while maintaining the battery state of charge (SOC) within limits.In [15], power divider laws which express branch power flows using the bus power injections in transmission networks are derived.In [13], a linearized power flow model based on piecewise linear approximation of the branch active power losses is proposed for solving the optimal power flow (OPF) problem in a radial distribution network.Though the proposed method is shown to have good convergence properties, the reactive power flow in the branches are neglected which makes it not suitable for voltage regulation and reactive power provision from flexible resources.One of the contributions of this paper is to extend the OPF method in [13] to include expressions for branch reactive power flows.The linear optimal power flow method developed in this paper is based on relaxation of branch active and reactive power losses and it is called branch losses relaxation (BLR) based LOPF.As compared to [13], the available voltage measurements or estimation values at each time step are utilized in BLR-LOPF to recompute the linear power flow matrices to achieve high accurate power flow solution and applicable for voltage regulation problems.
Though the flexibility from BESS and heat pumps only are investigated for demand response in this paper, the approach can be extended to other types of flexible resources such as electric vehicles.The proposed control method will be helpful to harness the available solar power and facilitate the customers with flexible resources in the LV network to participate in demand response.Over a receding time horizon, the proposed flexible resources control (FRC) computes power setpoints for all the flexible resources, taking into account key factors including (i) price signals from the electricity market [16], (ii) predictions of PV generation and loads, and (iii) the slow dynamics and/or capacity constraints of residential HP and BESS, and using network state estimates from an observability module (state estimator).The method can be used by the distribution system operator (DSO) for costeffective network power management.
The major contributions of this paper are 1) A new linear power flow method based on relaxation of branch power losses with improved accuracy in power flow calculations, and 2) A predictive control method for the control of flexible resources such as BESS and HP to provide demand response.
This paper is organized as follows.The architecture of the proposed predictive control method is discussed in Section II.In Section III, formulation of the models of linear power flow, BESS and HP are presented.The proposed MPC control method using LP technique is explained in Section IV.Simulation studies conducted on a Danish LV network using the proposed control are presented in Section V, and Section VI concludes the paper.

II. PREDICTIVE CONTROL OF FLEXIBLE RESOURCES
The major challenges of a DSO in the operation of MV and LV distribution networks with high DER penetration are voltage regulation and network congestion.To address these problems, a hierarchical control structure with three levels of control is proposed in this paper.The top level control consists of an MV network controller (MVNC) located at the HV/MV substation, the middle level control is exercised by aggregated LV network controllers (LVNC) at each MV/LV substation and the bottom level control (called as local control) is provided by individual device controllers.The main task of MVNC is to do the day-ahead planning of power dispatch at the MV network.Based on the forecasted generation and load, the MVNC estimates if there will be network congestions in the following day, computes the flexible power required from flexible resources at each MV bus and purchases the flexibility from day-ahead market.The MVNC then communicates the aggregated power setpoints to the LVNC [6].The proposed FRC which is part of the LVNC controls the flexible resources at the LV network such that the power setpoints from the MVNC are tracked, the LV network losses are minimized, the flexible resources are operated within their operational limits and the LV network constraints (network voltage and current limits) are satisfied [7], [9].The local controllers of individual flexible resource controls the power consumption of the equipment based on the received setpoints from the FRC [17].A typical distribution network with the above control architecture is shown in Fig. 1.In this paper, the focus is on the design of an predictive control algorithm for LVNC level with the layout described in the following subsection.

A. Layout of the Proposed Control Method
The layout of the proposed predictive control of flexible resources in a LV network is shown in Fig. 2. The operating time of each block can be identified by its color.It consists of four layers which are briefly explained as follows.
1) Physical layer: This layer consists of the LV network, flexible resources such as BESS and heat pumps, data hub and weather data resource center.The electricity consumption data of customers in Denmark are retrievable from data hub [18] for load forecasting and the weather information along with historical data are used for generation forecasting.The smart meters at customer locations measure parameters such as active and reactive powers, voltages and currents which are used by the distribution state estimation (DSE) algorithm to estimate the bus voltages.
2) Observability layer: This layer contains forecasting module, network state estimator (DSE) and observer of states of the flexible resources.Forecasting module provides the hourly predictions of the loads and PV generation aggregated at each LV bus.The proposed FRC depends on the state of the network which requires an accurate state estimation from network measurements.It is to be noted that state estimation in LV distribution networks could be possible with increasing deployment of smart meters at customer premises.As an example, more than 60% of Danish customers already have smart meters installed and by the year 2020, all Danish electricity end-users are obliged to install smart meters [18].A DSE algorithm which uses bus power injection measurements at critical buses and forecasts of load and generation at all other buses as pseudo-measurements could be employed to compute the bus voltage magnitudes [19].State observer block calculates the internal states of each flexible resource (state of charge of the BESS and floor temperature of the heating system) based on measurable quantities (voltage and current of the BESS and water and room temperature of the heating system) [17].
3) Control layer: The proposed FRC is part of the threelevel hierarchical control structure and it receives aggregated hourly power setpoints from the MVNC and hourly price signals from the electricity market.FRC computes the power setpoints of individual flexible resources by solving the proposed BLR-LOPF every 10 min.The local control comprises of the individual device controller, for example, the BESS charging/discharging unit which does the battery inverter control based on the power setpoints [6], [20].
4) Commercial layer: The day-ahead and regulation power markets (RPM) are involved in the proposed control structure to which the prosumers who own flexibile resources submit their bids through a commercial aggregator [20].
In this work, it is assumed that ideal forecasting, DSE, state observer and MVNC modules are in place and their design is not discussed.Instead, our primary focus is the control of flexible resources at the LV network which solves an optimization problem formulated in Section IV in a recedinghorizon fashion based on a linear power flow model of the network and linear models of the BESS and heating system which are formulated in the next section.

III. MODELING OF THE DISTRIBUTION NETWORK
In this section, the modeling of various elements of the network are explained.

A. Model of Distribution Network and AC Power Flow
In the modeling of the distribution network, matrices and vectors are represented using bold font; (•) T denotes matrix/vector transpose; |•| denotes the magnitude of a variable; (•) denotes a complex phasor; (•) and (•) denote the real and imaginary parts of a complex number respectively; j := √ −1 denote the complex unit; (•) * denotes complex conjugate and arg(x) denotes the phase of a complex number x.A balanced three-phase AC distribution network of radial topology is considered in this work.Such a network can be represented as an oriented acyclic graph G := (N 0 , L), specifically a root oriented tree.The vertices of the tree, denoted by the set N 0 := {0} ∪ N , represent the N + 1 network buses which consist of the set of N buses denoted as N := {1, . . ., N } and the slack bus (representing the secondary side of the LV distribution transformer) designated as 0. The edges of the tree denoted by the set L ⊂ N 0 ×N 0 represents the network branches with arbitrary numbering, which are oriented outwards from the slack bus [21].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.Let vn , īn ∈ C denote the phasors for the bus-to-ground voltage and current injected at each bus n ∈ N 0 whose magnitudes (|v n | , | īn |) are denoted simply as v n and i n respectively.As the slack bus voltage is taken as reference, arg(v 0 ) = 0 and v0 = v 0 .The complex bus-toground voltages and currents injected at the N non-slack buses are Let the vectors p ∈ R N and q ∈ R N contain the values of net active power and reactive power demand (load minus generation) at each bus n ∈ N .Let us define a vector the slack bus voltage magnitude v 0 .Let us define a diagonal matrix Vd = diag(v).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.The branch currents īb := [ īb1 , • • • , ībN ] ∈ C N ×1 can be computed from the bus currents ī as follows where, M r ∈ R N ×N is the reduced 1 bus injection to branch current matrix with binary values [22] and it is the inverse of the reduced bus incidence matrix (Note: The bus injections are considered negative).The branch voltage drops can be computed using Ohm's law from the bus current injections as vb = (R db +jX db ) īb .The bus voltage phasors v = v s −M T r vb can be expressed in terms of branch currents as It is to be noted that in (2), the matrix M T r is multiplied with individual branch voltage drops to get the corresponding voltage drops with respect to the slack bus.Using the matrix M r , the branch power flows at the sending end of each bus can be expressed as sb = M r (s + sl ), where, s ∈ C N is the vector of bus complex power injections, sl ∈ C N is the vector of complex branch power losses.From the branch currents and impedance, the complex power losses of the branches are 1 The column pertaining to the slack bus current is removed.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited.Content may change prior to final publication.and the complex branch power flows at the sending end of each bus are given by The equations ( 2)-( 4) can be written recursively in the form of DistFlow equations introduced in [23] for radial distribution networks, also called relaxed branch flow model in [24].The above equations exactly describe the nonlinear ac power flow problem and are difficult to use in an optimization framework without approximations [14].

B. State of the Art of Linear Power Flow Methods
The proposed control problem requires a linearized version of the original AC power flow equations.Linear power flow equations based on LinDistFlow method are presented in [12].The network losses are neglected in those power flow equations due to which the calculation of substation powers may be inaccurate.Another method for linear power flow is proposed in [8] which requires nominal voltage phasors, bus active and reactive powers for linearization.Though the above method is accurate for operating points close to nominal values, getting accurate measurements/estimates of voltage angles, bus active and reactive powers to re-linearize the system will be challenging.The method for piece-wise linear approximation of branch power losses is presented in [13].However, the reactive power flows are neglected, which may result in inaccurate calculation of reactive power injection from flexible from flexible resources.

C. Proposed Linear Power Flow Method
In the proposed method unlike [13], the reactive power flows in the branches are also taken into account in the formulation and the linear model is updated at each timestep based on the estimated bus voltages.The power flow linearization will be carried out at each step-time κ, using the latest voltage measurements/estimation.The aim is to express the bus voltages and branch power flows as linear functions of bus active and reactive powers (net power demand) through the following relations where, B vp , B vq and k v are constants at each time-step, F pl and F ql are linear functions of p and q which approximate and bound the actual branch active (p l ) and reactive (q l ) power losses respectively from above.The derivation of these expressions are done in the next subsections.In radial LV distribution networks, the R/X ratio is typically greater than 2.0 and the following assumptions are valid [25].
A1) The complex bus voltage phasors are approximately equal to the bus voltage magnitudes as the imaginary part of bus voltage phasors are negligible.A2) The linearization of power flow equations holds good for each time-step.1) Linear Approximation of Bus Voltages: Let us assume that voltage magnitudes of the N buses are measured at each sample time, denoted as 1) can be expressed in terms of the bus active and reactive powers as follows. īb where, the bus current injections are calculated as ī ≈ V −1 dm (p + jq) * .The real part of the expression in ( 2) is where, ∆v = v − vs .Using (8), above expression can be written in terms of bus active and reactive powers as below.Applying the first approximation stated above (i.e., v ≈ v), (10) becomes 11) with ( 5), the constants are defined as follows The above expression for bus voltage magnitudes is linear (as the matrices B vp and B vq are constants at each sample time) with respect to the bus active and reactive powers.
2) Piece-wise Linear Approximation of Branch Power Losses: Substituting (8) in (3) and expanding the terms, we obtain where, M rv = M r V −1 dm .In the approximation shown in ( 12), the active power losses p l ∈ R N ×1 + are separated into two parts (p lp and p lq ) based on the contribution of bus active and reactive power to the branch active power losses respectively such that p l = p lp + p lq .In a similar way, the branch reactive power losses are also expressed as q l = q lp + q lq .It is to be noted that the individual branch power loss approximations in ( 12) and ( 13) are quadratic with respect to either p or q.In order to linearize these expressions, we can construct piece-wise linear functions for a short range of branch current magnitude i b such that the approximated power losses bound the true losses from above by a small amount.In this way, a polyhedral relaxation of the quadratic bus powers [26] can be constructed as follows.Let the following approximation p lp1 ≈ B p1 p + d p1 be valid for the range i b ∈ [0, I b1 ] and p lp1 ≈ −B p1 p + d p1 for the range i b ∈ [−I b1 , 0] where, B p1 ∈ R N ×N and d p1 ∈ R N are constants at each sample time, such that p lp1 ∈ R N ×1 + is linear with respect to p. Similar terms can be derived for the whole range of i b as follows.Please refer to Appendix A for illustration of the above concept for a small distribution network.Let the set which are calculated as follows where, I b max is the maximum branch current magnitude.Let us define the following K diagonal matrices with branch resistances (B pk,k∈K ) and reactances (B qk,k∈K ) as follows.
where, 1 N denotes the N × 1 vector with all ones.Let us also define the column vectors given below.
From the above diagonal matrices and column vectors, let us define a linear function f pk (x) = B pk x + d pk .Using the above function, the approximate branch active power losses due to bus active power (p lp ) can be calculated from the below expression.
where, i ∈ N and e i denotes a unit vector in R N whose i th element is 1 and rest are 0. The approximate branch active power losses due to bus reactive powers (p lq ) are given by Similar expressions can be written for the branch reactive power losses using the piece-wise linear function f qk (x) = B qk x + d qk as given below.
The expressions in (17) define the function F pl in ( 6) and the equations in (18) define F ql in (7).
The above approximate expressions for branch power losses are expressed as relaxations (inequalities) in the linear OPF formulation to define an upper bound for the power losses and are provided in Section IV.Although the proposed linear power flow method is derived for balanced three-phase system, it is extensible to unbalanced distribution networks.In this case, there will be the three quantities for all the network variables (voltages, currents, active and reactive powers) to represent each phase separately [8].

D. Modeling of the Flexible Resources
In this subsection, the linear models of the BESS and heating systems which are used in the proposed linear OPF method for demand response application are explained.
1) Model of the BESS: In this work, a linear battery model developed in [27] which is a modified version of the kinematic battery model (KiBAM) proposed in [28] as given in ( 19) is adopted.The KiBAM consists of two charge capacity wells separated by a valve such that the charge is accessible from the first charge well.The battery model has two states (x bs = [q 1 q 2 ] T ) to model the available charge (q 1 ) and bound charge (q 2 ) respectively and the battery state of charge is SOC = q 1 + q 2 .The constraints on the battery charges are q 1 ≤ c w and q 2 ≤ 1 − c w .
The parameters c w and c r indicate the width of the first charge well and conductance of the valve with typical values of 0.1 pu and 0.001 pu respectively for Lithium-ion battery [27].The parameter C bs is the battery capacity in kWh.The charging and discharging efficiencies of the BESS are assumed equal in this work.The above model can capture the short-term battery dynamics [28] based on the discretization time-step and it is assumed that the power flow from/to the battery (p bs ) is  2) Model of the Heating System: In this work, we consider modeling of water-based HP which can provide high flexibility in power consumption [29].The HP transfers heat to the room through hot water flowing in pipes cast on the floor.A linear third-order model of a residential heating system with HP developed in [30] is used in this work, and its details can be found in [9].The continuous-time state-space representation of the heating system is provided below.
where, x hs = T r T f T w T are the room, floor and water ), p hs is the electrical power input to the HP in kW, T a is the ambient temperature in • C and ρ s is the solar irradiation.The typical values of the parameters of the heating system used in the simulations are provided in Table I.The time constant of this third-order system calculated from the dominant eigen value is about 7.5 hours.It is to be noted that zero-order hold method is used for discretizing the continuous-time models of BESS and heating system and therefore it is step-invariant [31].

IV. PROPOSED FRC BASED ON BLR-LOPF
In this section, the proposed optimal control method for flexible resources is derived.Let us assume that there are S number of BESS and H number of heat pumps connected to a LV network.Let the vector p g denote the PV generation at each node (the values are zeros at buses where PV is not connected) and the vectors p d and q d indicate the active and reactive power demand at each bus respectively.In this work, it is assumed that the PV inverters are operated at unity PF, hence q g = 0.The bus power injections p and q are calculated from the following.
where, C ss ∈ R N ×S and C hh ∈ R N ×H are binary matrices which map the BESS and heat pumps to the LV bus to which they are connected.The proposed optimal control is explained below.Let us define a variable κ = {1, • • • , N p } where, N p is the prediction horizon.
The power setpoints of the flexible resources is obtained by solving the optimization problem in (22).
The objective function (J) defined above has four components, which are the error in aggregated active power setpoint tracking (p se ), total active power loss of the network where, p l = p lp + p lq , sum of absolute value of reactive powers from all BESS (q bs,a ) and sum of auxiliary variables related to the branch reactive power losses (q tp and q tq ).The objective function is subjected to the following constraints.
The above constraints are used to realize the l 1 norm of (p κ sr − p κ s ) using the dummy variable p κ se .In the following constraints, the vectors p and q are given by the expressions in (21).The power balance equations are written as equality constraints provided below.
The bus voltages are constrained linearly to the bus active and reactive powers as provided below.
The following inequality constraints realizes the max operator in (17a) and (17b) for the branch active power loss components (p lp due to bus active powers and p lq due to bus reactive powers).
The branch reactive power loss components (q lp and q lq ) are also expressed as inequality constraints, but the realization of max operator in (18) requires auxiliary variables q tp and q tq .It is because, unlike the branch active power losses which are minimized via the objective function, the relaxed expressions of q lp and q lq are not "tight".Hence the auxiliary variables are introduced which will force the solver to tighten the upper bounds of q lp and q lq [32].The branch current and bus voltage limits are specified by the following box constraints where, B m = k m M r V −1 dm and k m is a scaling factor based on the branch power factor.
The state-space representation of the BESS are provided in the following equality constraints [27] x κ+1 bs = Φ bs x κ bs + Γ bs p κ bs (22l) The bounds of the active powers, reactive powers, SOC and internal states of the BESS are provided in the below constraints.p bs,min ≤ p κ bs ≤ p bs,max q bs,min ≤ q κ bs ≤ q bs,max x κ bs (1) ≤ c w (22r) The state-space representation of the heating system and the bounds of the heating system states are specified in the constraints below [9].
x κ+1 hs = Φ hs x κ hs + Γ hs p κ hs (22t) where, 0 N denotes the N × 1 vector with all zeros; c ps is the penalty cost for not tracking the aggregated power setpoint; c pl is the cost of LV network power loss; the vector c qbs contains the cost of reactive power provision from each BESS; the vectors q tp and q tq are auxiliary variables on constraints (22h) and (22i) to enforce an upper bound on branch reactive power loss q l ; the vector i b,max contains the branch current limits.The battery degradation is controlled by specifying limits for the SOC of the BESS in the range of 0.1 to 0.9 pu.The above optimization problem explicitly considers the bounds of room temperature using the constraint (22u) to ensure the thermal comfort of the customers and the network voltage limits in the constraint (22k) to ensure the voltage quality of supply.Hence, the flexible resources are managed in such a way that the quality of service is not compromised while harnessing the flexibility from the BESS and HP.
The weights of each component in (22a) varies at each steptime and are determined as follows.The weight c ps is set based on the price for regulating power imbalances [20], because it would be the price incurred to purchase balancing power from alternative sources.The weight c pl is chosen based on the hourly electricity price signals from the day-ahead market [33].The weight c qbs can be chosen based on the market price for reactive power, if applicable.In this work, c qbs is calculated based on the electricity price and inverter losses of BESS as discussed in [34].The procedure to solve the proposed predictive control method is stated below.
Algorithm 1 Proposed FRC based on BLR-LOPF 1: Get the state vector [v m , x bs0 , x hs0 ] from the observer blocks.2: Calculate B vp , B vq , B pk,k∈K , B qk,k∈K using v m .3: Solve the optimization problem defined in (22a) for the prediction horizon N p .4: From the optimal solution u, send only the power setpoints [p * bs , q * bs , p * hs ] for the next time-step to the local controllers.Go to step 1.
In this work, the PV inverters are considered to be operated at unity pf to show the worst-case scenario in network operation and emphasize the control of both active and reactive power of flexible resources such as BESS for maintaining the bus voltages within limits.In practice, however, the PV inverters may have to provide reactive power support based on grid code, PV inverter size and network conditions.Though it will reduce the active power capability of the PV inverters as their volt-ampere ratings are fixed, it will aid the proposed BLR-LOPF algorithm to maintain the network voltages within It is to be noted that when the states of a certain flexible resource are out of bounds, the flexibility available in that resource is zero, hence its power generation/consumption is not dispatchable by the DSO.For example, if the room temperature of a heating system participating in the proposed demand response scheme is above the comfort range, the HP is not operated even though its power consumption may be desired at that particular time.Only after the room temperature returns back to the comfort range, it can be utilized in the demand response scheme.If few flexible resource owners decide to opt out of the scheme for some reason, the overall flexibility of the proposed scheme may be reduced.However, the flexible resources will be able to rejoin the scheme later, if their states are within the desired range.In this work, the power levels of flexible resources are assumed to be continuous values dispatched by the DSO and the control process is assumed to be deterministic.The development of control strategies for DERs with uncertainty and discrete levels of output powers are left as a future research activity.
The proposed algorithm can be applied for unbalanced distribution networks in which the flexible resources may be connected as single-phase loads.This is possible by considering a set of three variables per phase and per LV bus to treat phase separately [8].The power setpoints of single-phase flexible resources will have to be computed on per phase basis.However, the detailed procedure and formulation of OPF for unbalanced networks is beyond the scope of this paper.

V. SIMULATION STUDIES
In this section, the simulation studies done for comparison with methods available in the state of the art and numerical validation of the proposed method formulated in Section IV are presented.Let us consider the 4-bus LV network shown in Fig. 11 with the following parameters.Let the impedance of first branch which represents the transformer impedance Z 12 = 0.0117 + j0.0327 pu and other branch impedances be 0.0102 + j0.0085 pu .Assuming the power injections at nodes 1 to 4 be −0.215−j0.09pu.The error in substation reactive power (Q s )  [12] and [13], while in the proposed method it is less than 0.5%.The error in substation voltage is about 3% [12], while in [13] and the proposed method it is less than 0.1%.The linear power flow method reported in [8] requires nominal voltage magnitudes and angles, bus active and reactive powers for linearizing the power flow equations, whereas the proposed method needs only the voltage magnitudes.Let us consider a change in active power injection from −0.2 − j0.09 pu to 0.1−j0.09pu.The change in error percentage in LPF [8] for P s and Q s are 23.2% and 32.0% respectively while for the proposed method the corresponding values are just 1.17% and 1.31% respectively.Hence the errors in estimation of power flows will be less in the proposed LPF compared to [8] when the network operating point (power injections at each bus) changes.In practice, it is difficult to measure/estimate values of voltage phase angles, bus active and reactive powers at all LV buses with good accuracy.Any errors in these values may have an impact on the results for the method proposed in [8].From the above results and discussion, it can be seen that the proposed LPF is more accurate in the estimation of network power flows.

B. Performance of Proposed BLR-LOPF
In order to assess the performance and accuracy of the proposed BLR-LOPF it is compared with the second-order cone relaxation programming (SOCP) based OPF [24] and the results are shown in Table III.Both the OPF methods are applied for solving the OPF of the network shown in Fig. 3 under typical network load conditions using Matlab on a desktop computer representative of current standards.It is to be noted that the SOCP based OPF is exact [24] and the results are global optimal points to the original AC OPF problem.The comparison results show that the proposed BLR-OPF is accurate though it slightly overestimates the branch power losses and underestimates the branch voltage drops.The execution time of proposed OPF is much less compared to the SOCP based OPF.

C. Simulation Studies of Proposed BLR-LOPF on Danish LV Network
A LV radial feeder with 17 nodes connected to residential customers shown in Fig. 3(a), which is a modified version of LV network at Askov, Denmark is used for this study.The ratings of the solar PV, BESS and heat pumps at each bus are provided in Fig. 3 Five simulation cases are studied to validate the significance of the proposed predictive control in demand response application and economic operation of the network.All the simulations are done with a time-step of 10 minutes for a period of 24 hours.In the below simulation cases, it is assumed that the device controllers (BESS and HP) which receive the power setpoints from LVNC are able to operate the respective devices without any additional constraints/limitations.The BESS at all nodes are assumed to have a minimum SOC of 0.2 at the beginning of the simulation.It is also assumed that the updated system states are available for the LVNC through proper communication channels at each time-step of 10 min.

(b). The plot of solar irradiance on a typical
0885-8950 (c) 2018 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

D. Case 1: Without Proposed Control of Flexible Resources during Summer
This simulation case is carried out to show the BESS operation without the proposed control method.In this simulation case, the BESS at each node is charged when the net power injection (difference between PV generation and load) is positive and the BESS is discharged if it is negative, provided the SOC is within the allowed limits.The heat pumps are not operated as the ambient temperature is high and there is no heating requirement.The plot showing the SOC of the BESS, bus voltages and transformer active power (p s ) is shown in Fig. 5.As seen from the figure, the BESS starts charging soon after solar PV start producing power resulting in non availability of flexible power at peak PV power generation.Due to unavailability of flexible power from BESS, the bus voltages at the end nodes are above the limit of 1.1 pu and the peak transformer loading is close to 100% during peak PV power.

E. Case 2: Flexibility from BESS during Summer
In this case, the proposed control method is applied to harness the flexibility of BESS as follows.The aggregated active power reference (p sr ) is set to limit the peak power infeed from solar PV to 75% of the transformer capacity.The objective here is to utilize the flexibility from the BESS during peak PV power and recycle the power during evening peak load period, thereby minimizing the net power exchange with MV network.The simulation results with a prediction horizon of 4 hours are shown in Fig. 6.The reference LV substation power is tracked with minimum error during peak PV power by postponing the charging of BESS to mid noon.All the BESS are discharged during evening peak demand minimizing the power consumption from the LV substation as the electricity prices are high during this period.The reactive power absorption from BESS at R17 can be seen in Fig. 6(d) to regulate the network voltages.As a result, all the bus voltages are maintained within the limit of 1.1 pu.

F. Case 3: Without Proposed Control during Winter
Simulation is done for a typical winter day to show the utilization of BESS for optimal operation of the network.The results of the simulation are shown in Fig. 7. Similar to Case 1, the BESS is charged/discharged based on net power injections at each bus and the HP are operated using a hysteresis control based on the room temperature bounds.The results of the simulation are shown in Fig. 7 and it can be seen that the bus voltages are slightly lower than the allowed limit of 0.9 pu.

G. Case 4: Flexibility from BESS during Winter
Simulation is done for a typical winter day to show the utilization of BESS for optimal operation of the network.Similar to Case 3, the the HP are operated using a hysteresis control based on the room temperature bounds and the BESS is operated based on the proposed FRC.The results of the simulation are shown in Fig. 8 and it can be seen that BESS are utilized well compared to Case 3 to meet the network control objectives.However, similar to case 3, the network voltages are slightly lower than the limit of 0.9 pu.

H. Case 5: Flexibility from BESS and HP during Winter
Simulation is done for a typical winter day to show the utilization of both BESS and HP for optimal operation of the network.The room temperature limits for all the heating systems are set in the range 19 to 23 • C. The simulation results are shown in Fig. 9.It can be seen in Fig. 9 that the BESS and HP are operated to track the aggregated LV power setpoints within the constraints of SOC and the comfort range of room temperatures set by the user.The HP are operated such that the room temperatures are close to the upper limit (23  to the lower limit (19 • C) when the prices are high.Also, the network voltages are maintained within the limit of 0.9 pu as seen in Fig. 9 (b).

I. Calculation of Flexibility from BESS and HP
As discussed in [35], the flexibility of resources can be assessed in two ways, the first one based on comparing the power consumption at each time instant without and with flexibility control and the second method based on the costs incurred without and with activation of flexibility.Let us define flexibility index (φ) of each flexible device [35] over an operating period of N time-steps to be where, p f i,wc and p f i,c are the power consumptions of i th flexible resource without and with control respectively, p f i max is the maximum power consumption.The values of p f will be negative for BESS during discharge.The values of φ are in the range [0, 1].The computed flexibility indices φ bs and φ hs of the BESS and heating systems respectively over a period of 24 hours are shown in Fig. 10.As seen from Fig. 10, flexibility obtained from BESS is slightly higher in summer due to high PV generation compared to winter.The economic savings index (ESI) [35] due to utilization of flexibility can be computed from the equation below.
where, J wc and J c are the performance indices (total costs calculated from (22a) and summed up for 24 hours) of the network without and with proposed control.As seen from Table IV, the cost savings are more in summer as the solar power can be effectively stored by utilizing the flexibility from BESS and then it can be discharged at evening peak load times.

J. Summary of Simulation Studies
From the simulation results of Case 2, it can be seen that the proposed control method can provide demand response of BESS to reduce the peak transformer loading by 25% on a typical summer day.From Case 5, it can be seen that proposed control can enable optimal operation of BESS and HP to track the reference aggregated power while satisfying the customer comforts and device constraints during winter.In both Case 2 and 5, the results show that the BESS and HP are operated to consume power during off peak hours when the electricity prices are cheap and discharge of power from BESS/less power consumption by HP when the electricity prices are high resulting in economic operation of the LV network.

VI. CONCLUSIONS
In this paper, a predictive control method is proposed for the aggregated low-voltage network control to effectively utilize the flexible resources such as BESS and HP for demand response application.The results of the simulation studies show that the proposed method can find setpoints for the flexible resources for economic operation of the network based on the electricity prices, forecasts of generation and load and system states.As a future work, a distributed control method to solve the above optimization problem with minimum exchange of system information among the controllers, thereby ensuring data privacy of the flexible resource owners will be pursued.

APPENDIX A APPROXIMATION OF BRANCH POWER LOSSES
The notations described in Section III is shown for a simple 4-bus LV network in Fig. 11.It is to be noted that the labeling of the branches is arbitrary.The relation between branch power flows and bus powers can be written as given below.In a similar fashion, p lq2 , q lq2 and q lq2 can be calculated.For each branch, four expressions of branch power losses such that of (26) have to be evaluated.The max operator in these expressions are converted to inequality constraints in linear OPF resulting in 8k × N constraints where, k is the number of branch current partitions and N is the number of branches.

APPENDIX B PARAMETERS OF THE LV NETWORK
The cable parameters of the LV network shown in Fig. 3 (a) are provided in Table V.The resistance and reactance values are per phase quantities.The total distance from bus R1 to R7 is 0.69 km and the distance from R3 to R17 is 0.45 km and from R6 to R17 is about 0.3 km.

Fig. 3 .
Fig. 3. LV network considered for studies (a) schematic, (b) ratings of the PV and flexible resources

Fig. 4 .
Fig. 4. Plot of profiles during winter and summer: (a) solar irradiance measured on a clear sky day, (b) typical ambient temperature (c) hourly load profile of a residential house and (d) typical hourly price signals of regulating power market (RPM) for down regulation and day-ahead market (DAM).

Fig. 5 .
Fig. 5. Plot of case 1: without proposed control during summer (a) SOC of the BESS (b) bus voltages at the end nodes of the network (c) aggregated LV power flow.

Fig. 6 .
Fig. 6.Plot of case 2: flexibility from BESS during summer (a) SOC of the BESS, (b) bus voltages, (c) aggregated power tracking and (d) active power setpoints to BESS at R17.

Fig. 8 .
Fig. 8. Plot of case 4: flexibility from BESS during winter (a) SOC of the BESS, (b) bus voltages, (c) aggregated power tracking, (d) active power setpoints to BESS at R17 and HP at R16 and (e) room temperatures of residential heating systems.

Fig. 9 .
Fig. 9. Plot of case 5: flexibility from BESS and HP during winter (a) SOC of the BESS, (b) bus voltages, (c) aggregated power tracking, (d) active power setpoints to BESS at R17 and HP at R16 and (e) room temperatures of residential heating systems.

Fig. 10 .
Fig. 10.Plot of flexibility indices of (a) BESS and (b) heating systems
This article has been accepted for publication in a future issue of this journal, but has not been fully edited.Content may change prior to final publication.Citation information: DOI 10.1109/TPWRS.2019.2898425,IEEE Transactions on Power Systems 0885-8950 (c) 2018 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited.Content may change prior to final publication.Citation information: DOI 10.1109/TPWRS.2019.2898425,IEEE Transactions on Power Systems 0885-8950 (c) 2018 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
In this work, the optimization problem is solved every 10 min, hence the minute variations in battery states are not captured by the discretized model.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited.Content may change prior to final publication.Citation information: DOI 10.1109/TPWRS.2019.2898425,IEEE Transactions on Power Systems 0885-8950 (c) 2018 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited.Content may change prior to final publication.Citation information: DOI 10.1109/TPWRS.2019.2898425,IEEE Transactions on Power Systems SUBMITTED TO IEEE TRANSACTIONS ON POWER SYSTEMS.THIS VERSION FEBRUARY 7, 2019.8 A. Comparison of Proposed Linear Power Flow with State of the Art 0885-8950 (c) 2018 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited.Content may change prior to final publication.Citation information: DOI 10.1109/TPWRS.2019.2898425,IEEE Transactions on Power Systems The calculated ESI values for different simulation cases are provided in TableIV.0885-8950(c) 2018 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.