Data-Driven Risk Analysis for Probabilistic Three-Phase Grid-Supportive Demand Side Management

Along with the emerging development of demand side management applications, it is still a challenge to exploit flexibility realistically to resolve or prevent specific geographical network issues due to limited situational awareness of the (unbalanced low-voltage) network as well as complex time dependent constraints. To overcome these problems, this paper presents a time-horizon three-phase grid-supportive demand side management methodology for low voltage networks by using a universal interface that is established between the demand side management application and the monitoring and network analysis tools of the network operator. Using time-horizon predictions of the system states that the probability of operational limit violations is identified. Since this analysis is computationally intensive, a data driven approach is adopted by using machine learning. Time-horizon flexibility is procured, which effectively prevents operation limit violation from occurring independent of the objective that the demand side management application has. A practical example featuring fair power sharing demonstrates the effectiveness of the presented method for resolving over-voltages and under-voltages. This is followed by conclusions and recommendations for future work.


Introduction
The fast growing share of distributed renewable energy sources (DRES) like solar photovoltaic (PV) and highly energy-intensive appliances and distributed energy resources (DER) such as heat pumps and electric vehicles (EVs) results in increasing uncertainties in the power flows in distribution networks, which challenges the distribution system operator (DSO) to keep the network operated within safe and secure operation limits [1][2][3][4].Due to the high uncertainties, a continuation of the current paradigm of infrastructure over dimensioning is expected to result in high future investment costs.To cope with this, DSOs can deploy control algorithms to resolve operational limit violations, e.g., using on-load tap-changers (OLTC) [5] or reactive power control of inverters [6,7].However, the application of such direct control actions can be highly expensive or otherwise ineffective.The consequence is that specific operation limit violations stay unresolved if the local controllers cannot effectively resolve the operation limit violation, i.e., there is a remaining operation limit violation (OLV).Alternatively, the DSO can invoke flexibility in the demand and supply of DER owners instead by using demand side management (DSM) applications [8].Due to high R/X ratio's typical in distribution networks, active power from DSM will have a more significant influence on the system states than a reactive power control [9].This flexibility can result from DER capable of advancing or deferring their starting time or changing the active power consumption or production during some time interval [10].Examples include various types of controllable appliances like time shiftable appliances (e.g., freezers or washing machines) and buffering appliances involving some form of storage such as batteries of EV and heat pumps [11].
DSM applications are expected to be operated by a market actor, i.e., aggregators or energy suppliers with limited to no insight in the network operation state.As such, they tend to be unaware about specific grid related issues and physical and geographic aspects of the network during operation.Although some are designed to resolve operation limit violations [12,13], their optimization objective is often not focussed on resolving specific operation limit violations occurring at a certain geographical location.This is especially true for low-voltage (LV) networks where uncertainty is even higher due to the non-aggregated load profiles [14].
These shortcomings can be solved by using a universal interface between DSOs and DSM applications, facilitating the exchange of information on specific (predicted) network issues, and ending user flexibility.In this case, state prediction of the network system states will be of high value to account for the probabilistic system states that will occur in the near future.For coping with the high stochasticity of DER and end user behaviour, a probabilistic approach is required to determine potential network risks [15][16][17].On the down side, probabilistic approaches often require some form of probabilistic power flow (PPF) calculations to evaluate the probability of having operation limit violations at certain locations in the network.Despite abundant attempts (e.g., [18][19][20][21][22]) to lower the computational complexity for PPF calculations for a large three-phase network, they still require a significant amount of time to complete.Performing PPF calculations for each 15-min time interval within a 24-h day-ahead period for all the networks it operates forms an enormous computational burden.Quite often, however, we are not interested in the full probability density functions produced by the PPF, but solely in the probability that a certain operation limit will be violated.Rephrasing along with control and management decisions are made based only on specific information of the operational limit violations (e.g., the probability of it exceeding a certain predefined threshold).As such, only a limited amount of information is required, which creates possibilities to speed up the computation time.
To this extent, in this work, we use machine learning to come up with predictions on whether the probability for operation limit violations exceeds a certain threshold and how DSM can bring the probability back to acceptable levels.Various machine learning techniques such as neural networks (NN) have attracted greater attention for applications in DSM.Often, their application can be found in load forecasting [23,24].Furthermore, various types of home energy management systems implement machine learning for decision making.For example, in Reference [25], a NN is applied for the scheduling of PV panels and a storage system.This work presents the application of NN to specifically deal with geographically dependent operational limit violations in distribution networks, which results in a probabilistic approach for time-horizon DSM using a universal applicable interface between DSO and DSM applications.Specifically, the probabilistic NN based analysis will indicate whether OLV are expected to occur at certain geographical locations in the network and with what probability.If this probability exceeds a certain threshold, the DSO will request the DSM application for flexibility to reduce the probability of the violation of operation limits at acceptable levels.Based on the probability of the system states of the three-phase nodes over a day-ahead (DA) or intra-day (ID) time interval, the three-phase sensitivity of the operational limit violation faced by the DSO with respect to the active power injection by customers is derived.The DSO will use this information to specify the expected operation limit violation and how its probability can be reduced towards the DSM application by using the universal interface.The main contributions of the paper are:

•
Time-horizon analysis of operation limit violations using a probabilistic NN based approach: if operational limit violations are expected with a certain probability, grid-supportive DSM is triggered by using a universal interface to reduce the probability back to acceptable levels; • Specification of a three-phase unbalanced network operation limit violations occurring with a certain probability over time and their sensitivity with respect to active power, according to the network model in Reference [26].This allows DSM applications to resolve network issues by optimizing time dependent flexibility, which is independent of the objective of the DSM application; • A demonstration will be given of a time-horizon DSM optimization resolving under and over voltages featuring fair power sharing among the different prosumers.
The remainder of the paper is organized as follows.In Section 2, the universal interface between the DSO and DSM applications is described, which allows for a generalized procurement of flexibility.Section 3 follows up with how the required amount of flexibility is triggered based on predicted system states and the network sensitivity, which is optimized according to the objective of the DSM application.Section 4 presents an alternative method based on a NN.Section 5 will present the overall DSM optimization while Section 6 gives a proof of concept using numerical simulations, which shows the applicability of the proposed approach on the thee-phase IEEE low voltage EU network.Lastly, Section 7 will elaborate further recommendations and conclusions.

Universal Procurement of Flexibility
This section briefly introduces the universal interface adopted in this work for establishing DSM using a probabilistic, time-horizon, grid-supportive methodology.A more in-depth discussion of this framework can be found in Reference [27].One of the main features of this interface is that it should work efficiently with all possible DSM applications in the field.In this case, the starting point of the interface is that the exchanged information should be non-iterative and integrable in the widest selection of optimization algorithms.Therefore, for each time interval, the DSO can specify what flexibility it requires depending on the geographical locations of the flexibility within the network.After receiving this information, the DSM application can allocate the required flexibility depending on its own optimization strategy and send the results back to the DSO.The specified information by the DSO on the required flexibility consists of lower limits for certain linear combinations of changes in active power injection at the geographical locations of the customer connection points and their corresponding linear gains.From this, linear constraints can be constructed, which can be considered by virtually any DSM application.These gains are formed by the sensitivity matrix of the network, which specifies the linearized sensitivity of the system states with respect to changes in active power, as detailed in Section 3.2.As such, independent of what objective a DSM application pursues, it will be able to resolve the OLV by including constraints on the before mentioned sensitivity of active power injection towards the OLV the DSO faces.This way, the DSM application is enabled to consider the operation limits at each geographical location in the distribution network in its optimization.Depending on the DSM application in place, additional information can be included such as for bidding or pricing information.In the proposed framework of Reference [27], the procurement of flexibility takes place in two stages: (1) time-horizon flexibility based on both probabilistic power flow and machine learning during the DA/ID preventive planning phase, which is the main topic of this paper and (2) real-time flexibility based on state estimation during the corrective execution phase, as described in Reference [27].
Preventive time-horizon flexibility will be procured in on a DA/ID basis.Preventive DSM based on predictive load forecasting [28][29][30] will be of high importance to account for the system states and possible operational limit violations that might occur in the near future.One of the motivations for this is that prosumer flexibility in energy production or consumption is expected to involve complex time-dependent interdependencies.This means that flexibility provided within some time interval might need to be compensated by flexibility in another time interval.Furthermore, to prevent excessive life time degradation, DSOs will prefer to not switch certain controllers very frequently.This can, for example, hold true for OLTC tap position adjustments.As such, an optimization over a certain time horizon of the network system states and the available flexibility is of high importance to resolve OLVs that are expected over time.For this, probability density functions (PDF) for each of the system states can be obtained for the time interval of the DA/ID optimisation period from where the probabilistic Energies 2018, 11, 2514 4 of 18 system states can be obtained by using probabilistic power flow calculations.DSOs can make important decisions with this information for adjustments of set points of the previously mentioned local control capabilities.Whenever there are any OLV after adjustments of the local controller setpoints, the DSO can trigger time-horizon flexibility from DSM applications by specifying the specific OLV at hand along with the sensitivity regarding active power injection at each geographical network location.
Within the corrective execution phase, real-time DSM is based on actual state estimation (SE) of the network.It is expected that SE capabilities must be expanded from transmission to distribution networks to gain insight in the network system states and allow for the control applications, which is mentioned earlier [31].Within the real-time execution phase, real-time insight in the actual system states will be crucial to trigger these control applications as well as real-time flexibility from DSM applications, which is described in Reference [27].

Probabilistic Grid-Supportive Flexibility
As stated, the DSO can invoke flexibility in active power injection from DSM applications to address OLVs.Based on the real-time and predictive monitoring capabilities, the interface between the DSO and the DSM application will facilitate a generalized information exchange on the OLV the DSO is facing with the sensitivity for changes in active power at specific geographical locations.As discussed, a newly proposed benchmarking method is proposed based on PPF calculations as well as a novel method using a NN.A flow diagram of the full functioning of the proposed methods covered in Sections 3-5 is shown in Figure 1.This figure illustrates both the benchmark approach and the NN approach.
Energies 2018, 11, x FOR PEER REVIEW 4 of 17 specifying the specific OLV at hand along with the sensitivity regarding active power injection at each geographical network location.Within the corrective execution phase, real-time DSM is based on actual state estimation (SE) of the network.It is expected that SE capabilities must be expanded from transmission to distribution networks to gain insight in the network system states and allow for the control applications, which is mentioned earlier [31].Within the real-time execution phase, real-time insight in the actual system states will be crucial to trigger these control applications as well as real-time flexibility from DSM applications, which is described in Reference [27].

Probabilistic Grid-Supportive Flexibility
As stated, the DSO can invoke flexibility in active power injection from DSM applications to address OLVs.Based on the real-time and predictive monitoring capabilities, the interface between the DSO and the DSM application will facilitate a generalized information exchange on the OLV the DSO is facing with the sensitivity for changes in active power at specific geographical locations.As discussed, a newly proposed benchmarking method is proposed based on PPF calculations as well as a novel method using a NN.A flow diagram of the full functioning of the proposed methods covered in Sections 3-5 is shown in Figure 1.This figure illustrates both the benchmark approach and the NN approach.

Probabilistic Prediction of Operational Limit Violations
In this paper, we assume an unbalanced radial network [26] consisting of N nodes including the slack node where there are M households tapping of the feeders offering flexibility in their active power consumption or production.The slack node is located at the root of the radial feeder and is assumed to be constant.To determine the required amount of flexibility in active power consumption or production at each node in the network within the preventive planning phase, the DSO performs a probabilistic analysis of the network system states for each time interval t of the DA/ID optimization.This can be either based on the before mentioned PPF or the NN approach presented in this paper.For benchmarking purposes, in this section, first the PPF approach will be discussed, which is followed by the NN based approach in Section 4.
As an input, the PPF calculation takes the power injection PDF of each household, which is based on predictions.The PPF calculation results in a PDF for the all the system states of the network.For example, the probability density of the voltage magnitude V t n,p occurring at node n, phase p, and time interval t is given by P V t n,p .Now suppose that the network has a lower limit and upper limit V l and V u for the voltage magnitude, then the probability of an undervoltage occurring at node n and phase p during time interval t can be calculated, according to the cumulative distribution function (CDF) F V t n,p .This is shown below.
Similarly, for the probability of over voltages, we can state that: Similar expressions can obviously be made for overloading or violation of voltage imbalances throughout the phases.However, in this paper, we focus only on under/over-voltages.Now, if the probability of having an under/overvoltage as expressed in Equations ( 1) and ( 2) exceeds a certain threshold value f .The DSO will opt to procure flexibility from the DSM application in order to shift the PDF of the voltage magnitudes such that the probability of having an under/overvoltage is brought down to acceptable levels.To formulate this as an easier condition to deal with, we first introduce the variables V t,l n,p and V t,u n,p .In this case, V t,l n,p is the voltage magnitude for which the CDF F V t n,p equals the acceptable probability threshold f .
In other words, there is a probability of exactly f that the voltage at node n and phase p at time interval t will be lower or equal to V t,l n,p .Similarly, V t,u n,p , can be expressed as: Rephrasing in other words, there is a probability of exactly f that the voltage at node n and phase p at time interval t will be higher or equal to V t,u n,p .In case of a discontinuous CDF, V t,l n,p and V t,u n,p take the next available value closest to the nominal voltage.We define the vectors V t,l p and V t,u p of all voltages V t,l n,p and V t,u n,p of all nodes n (apart from the slack node) in phase p and time interval t, according to the equations below. respectively, Energies 2018, 11, 2514 Now, the DSO will procure flexibility for under-voltages if any of the voltage magnitudes V t,l n,p are lower than the lower limit V l and flexibility for overvoltages if any of the voltages V t,u n,p are higher than the upper limit V u .After all, in both cases, the probability of having a voltage limit violation is higher than the probability threshold f .This can be expressed as the conditions below. and If one of these conditions holds, the DSO needs flexibility to reduce the probability at least to the acceptable probability threshold f , which is illustrated in Figure 2 and detailed in the next sections.
Energies 2018, 11, x FOR PEER REVIEW 6 of 17 higher than the upper limit  .After all, in both cases, the probability of having a voltage limit violation is higher than the probability threshold .This can be expressed as the conditions below.
, , ≤ and If one of these conditions holds, the DSO needs flexibility to reduce the probability at least to the acceptable probability threshold , which is illustrated in Figure 2 and detailed in the next sections.

Network Sensitivity Operation Point
To procure the right amount of flexibility from the DSM application, the DSO needs to know how much flexibility will be required to solve the OLV depending on the geographical location at which the flexibility is delivered.As stated in the introduction section, this work uses the linearized sensitivity of the network system states, which has been used effectively in other studies [27,32].The sensitivity of the probabilistic OLV informs the DSM application what linear combinations in active power flexibility can address the OLV.To obtain the linearized sensitivity of OLVs with respect to changes in active power injection ∆ , at node , phase , and time , the Jacobian matrix of the partial derivatives of the power injections  , with respect to the violated system states are derived.As an example, concerning violation of the voltage magnitude, the partial derivatives can be expressed by using well-known power flow equations as well as expressing the nodal power injection in terms of the nodal voltage magnitudes and network parameters (which are considered static and, therefore, left out of the notation).
It should be noted that if  < , the nodes at which no active power flexibility is available can be left out of the Jacobian.
From this point, the operation point of the partial derivatives is based on the outcome of the PPF calculation.The PPF calculation is in practice often and is completed by using Monte Carlo simulations, which results in discrete samples of the PDFs of the system states.In this paper, a pragmatic approach is taken to obtain an effective operation point from the discrete samples of the PDFs of the system states.Suppose that the Monte Carlo simulations for each time interval  result in a set  of discrete samples of the PDFs of the system states.For each sample  ∈  at time interval , we can compose the vector of all system states for each node  in phase .

Network Sensitivity Operation Point
To procure the right amount of flexibility from the DSM application, the DSO needs to know how much flexibility will be required to solve the OLV depending on the geographical location at which the flexibility is delivered.As stated in the introduction section, this work uses the linearized sensitivity of the network system states, which has been used effectively in other studies [27,32].The sensitivity of the probabilistic OLV informs the DSM application what linear combinations in active power flexibility can address the OLV.To obtain the linearized sensitivity of OLVs with respect to changes in active power injection ∆P t n,p at node n, phase p, and time t, the Jacobian matrix of the partial derivatives of the power injections P t n,p with respect to the violated system states are derived.As an example, concerning violation of the voltage magnitude, the partial derivatives can be expressed by using well-known power flow equations as well as expressing the nodal power injection in terms of the nodal voltage magnitudes and network parameters (which are considered static and, therefore, left out of the notation).
It should be noted that if M < N, the nodes at which no active power flexibility is available can be left out of the Jacobian.
From this point, the operation point of the partial derivatives is based on the outcome of the PPF calculation.The PPF calculation is in practice often and is completed by using Monte Carlo simulations, which results in discrete samples of the PDFs of the system states.In this paper, a pragmatic approach is taken to obtain an effective operation point from the discrete samples of the PDFs of the system states.Suppose that the Monte Carlo simulations for each time interval t result in a set K t p of discrete samples of the PDFs of the system states.For each sample k ∈ K t p at time interval t, we can compose the vector of all system states for each node n in phase p.
We define the subset K t,l p ⊂ K t p to be the set of samples k, which yields that the minimum value of V t,k p has a probability lower or equal to f , i.e., the minimum value of V t,k p is smaller or equal than the minimum value of V t,l p .minV t,k p ≤ minV t,l p (11) Similarly, the subset K t,u p ⊂ K t p is the set of samples k, which yields that the maximum value of V t,k p is larger or equal than the maximum value of V t,u p : Lastly, we define δ t,l p and δ t,u p as the average voltage angle vector, which is averaged elementwise over all complex voltages corresponding to all k being an element of the set K t,l , respectively, K t,u .Now, the operating point of the network for the Jacobian is chosen to be V t,l p ∠δ t,l p or V t,u p ∠δ t,u p , i.e., the voltage vectors for which the magnitude corresponds with the probability f and the angle is averaged over the phasors of the sets K t,l end K t,u .Note that the sensitivity will be different for under and over voltages, which are denoted as S t V,l and S t V,u .Averaging of the angles δ t,l p and δ t,u p is required since there is not a single angle that corresponds with the system states for the probability f .For the voltage magnitudes, there exists a unique relation, which is expressed in Equations ( 3) and ( 4).For the angles, such a unique relation does not exist.Therefore, we average over all the values in the sets K t,u p and K t,l p .Note that this is only used for determining a suitable operation point for the Jacobian and in no way aims to change the power factor of any appliance.
In relation to this, an important note should be made concerning the Jacobian.The partial derivatives forming the entries of the Jacobian differ depending on the relation between the reactive and the active power at each node and phase.This relation is determined by the appliance associated to the power injections.As an example, the power factor can remain constant for any change in the active power injection (therefore, changing the reactive power injection) or the reactive power can remain constant independent of the active power injection.This results in different partial derivatives where highly non-linear or discontinuous relations between active and reactive power will occur.This might complicate the process of deriving a suitable Jacobian.In the simulation results of this work, we assume the reactive power to remain constant.Lastly, the sensitivity with respect to the OLV is obtained by inverting the Jacobian, according to the equation below.

Constraints for Demand Side Management
To reduce the probability of a specific operation limit violation to acceptable levels, the DSO will specify the required minimum required change in any network system state.As flexibility for resolving operation limit violations most likely will result in shifting power consumption to another point in time, we also need to specify the available 'capacity' in time intervals where no operation limit violations are expected.If the OLV concerns nodal voltage magnitude violations, the DSO will specify the vectors ∆V t l and ∆V t u , which indicates the voltage with which the system states are exceeded and what capacity is available.We rearrange the elements of ∆V t l and ∆V t u , in general denoted as ∆V t , such that the vectors contain the elements for all nodes n (excluding the slack node) and all phases p ∈ {a, b, c}.
The individual elements of these vectors are given by the equation below. and In this case, V min and V max are the minimum and maximum voltage limits while V t n,p is the voltage before DSM at node n and phase p.
As a final step, the change in active power at time t for a certain node n, phase p, and time t is represented by ∆P t n,p where the elements for all nodes and phases together compose the vector ∆P t .
Currently, the sensitivity with respect to the OLV is obtained by inverting the Jacobian according to S = J −1 .herein this case, the linearized change in the network system states is calculated by multiplying the change in active power ∆P t n,p injection with the sensitivity matrix.For an OLV of the nodal voltage magnitude, the resulting change in the system states ∆V t n,p yields the following equation.

Neural Network-Based Grid-Supportive Demand Side Management
As discussed, performing the PPF is computationally costly despite the many studies in literature improving its efficiency.By carrying out DA PPF calculations for all time intervals in the DA period and all the networks it is managing, this might require significant computational power.Therefore, this work introduces an NN based approach to accurately approximate the findings derived in the previous section and, therefore, drastically speeding up the required computation times.After all, the interest in this case is not on the probabilistic system states but rather the need for flexibility.For this purpose, a multi-layer NN (NN) is introduced that only needs to be trained once for a particular distribution network.After training, evaluating the NN is considerably faster than performing the PPF, which results in a significantly reduced computational effort.The next subsections will respectively discuss the NN architecture and the training of the NN.

Neural Network Architecture
The NN is designed to prevent the costly PPF within the risk analysis of the DSO.To this extent, it needs to approximate the information on the operational limit violation at hand and its sensitivity with respect to changes in active power, which is detailed in Section 3.For this purpose, in this work, a regressive NN is applied to provide the DSM application with the required information by replacing the costly PPF.The overall architecture of the network is displayed in Figure 3.
For the PPF method described in Section 3, the inputs are the PDFs of each of the households while the outputs are formed by the required change in systems states ∆V t and the sensitivity matrix S t V .However, the full PDFs of the households are not very suitable to be used as an input for a NN since it normally expects a numeric input value rather than a continuous function specification.On a similar note, the output sensitivity matrix S t V has a number of entries equal to the square of the number of nodes, which makes it too large to be effectively approximated by a regression based NN.As a final point of concern, the required change in system states ∆V t is a non-linear function with a discontinuous derivative and, therefore, is also not very suitable for regression analysis.

Neural Network Architecture
The NN is designed to prevent the costly PPF within the risk analysis of the DSO.To this extent, it needs to approximate the information on the operational limit violation at hand and its sensitivity with respect to changes in active power, which is detailed in Section 3.For this purpose, in this work, a regressive NN is applied to provide the DSM application with the required information by replacing the costly PPF.The overall architecture of the network is displayed in Figure 3.For the PPF method described in Section 3, the inputs are the PDFs of each of the households while the outputs are formed by the required change in systems states ∆ and the sensitivity matrix To overcome these design challenges for the inputs of the NN, the PDF of the loading of each household can be fitted on a suitable well-known default PDF like a Gaussian, beta, or Weibull distribution.Specific features or shape parameters can be extracted like the mode, mean, median, and variance or the α and β or k and λ parameters for the beta and Weibull distributions.These shape parameters can be used as inputs for the NN.For the sensitivity matrix, although large matrix inversion is involved, determining the sensitivity based on a set of operation points of the network is a simple task and relatively computationally efficient.In addition, determining the required change in the system states ∆V t from the operation points is highly straightforward.Therefore, the NN is trained to not output the required change in the system states and corresponding sensitivity but rather approximate the operation point used for calculating the network sensitivity.
Mathematically expressed, these are V t,l p ∠δ t,l p or V t,u p ∠δ t,u p as introduced in Section 3, i.e., the voltage vectors for which the magnitude corresponds with the probability f and the angle is averaged over the phasors of the sets K t,l end K t,u .The network sensitivity S t V can be calculated based on Equations ( 9) and (13).Similarly, ∆V t is easily derived using Equations ( 15) and (16).
The simulation results presented in Section 6 of this work use data for the household loading PDFs that is suitable to be fitted on a Gaussian distribution.The mean µ t m and variance σ t m of the PDFs for each household m at time t are defined as the input variables of the NN.This means that the number of inputs for the NN will be equal to 2M.Similarly, the output variables together V t,l p ∠δ t,l p and V t,u p ∠δ t,u p will be of size 12(N − 1), (three phase nodes and four variables per node excluding the slack node).However, since the numerical values of the voltage magnitudes and angles are considerably different, it is better to split the NN in two separate networks of size 6(N − 1) since, this way, performance indicators such as the mean squared error (MSE) are used more meaningful.Lastly, in the average distribution network, we do not expect differences in voltage angles up to π radians.Therefore, it is advisable to shift the voltage angles with π radians to eliminate the transition between 0 and 2π.
As an alternative, one could consider working in Cartesian coordinates rather than polar coordinates.As a final step, experimental simulations will be required to determine a suitable number of hidden layers and neurons.For the experimental results presented in Section 6, it turns out that two or up to three hidden layers strike a reasonable balance between training time and accuracy by estimating the non-linear relation between the input and output variables accurately, which can be seen from the presented results in Section 6.The number of hidden neurons is highly dependent on the expected correlation of the input data and, therefore, should be determined experimentally for each network.

Neural Network Training
To train the NN with the architecture as described above, a large amount of historical data is required, which comprises the PDF of household loadings for many different situations.Based on this, the mean µ t m and variance σ t m can be derived together with the corresponding V as described in Section 3, which forms the training set for the NN.Clearly, it is important that the training set contains a sufficiently diverse number of situations that might occur in the network to make the system robust for unexpected events.Since the proposed NN architecture will have a significant size for larger distribution networks, it is advisable to perform the training of the NN on a graphics card to exploit the possibilities of parallelism.When doing so, from experimental results, backpropagation training using gradient derivatives and the steepest descent turns out to strike a good balance between performance and training time for the NN architecture proposed in this study.However, the convergence is very sensitive to the learning rate and, therefore, in this study, the gradient descent with an adaptive learning rate backpropagation is used, which is implemented in MATLAB (MathWorks Inc., 2018b (prerelease), Natick, MA, USA) [33].

Overall Demand Side Management Optimization
After the reception of the specification on the OLV and corresponding sensitivity from the DSO (either using the PPF or NN approach), the DSM application deploys an overall optimization of the available flexibility in active power injection depending on its optimization objectives and is constrained by the specified OLVs.As stated in the introduction section, the objectives can take many forms like local supply and demand matching [11], fair power sharing [34], or a market mechanism where global welfare is optimized [35].We define the set F of size M of all pairs of nodes and phases {m, q} at which flexibility is offered by prosumers.Independently of the optimization objective the DSM pursues, the general optimization function in Equation ( 19) can be defined alongside the corresponding constraints for voltage OLVs Equation (20).
subjected to: V,u ∆P t other constraints of flexible appliances (20) Examples of other constraints of flexible appliances can be found in Reference [36].Simultaneous over voltages and under voltages at different nodes of the network are (though unlikely because of the huge voltage difference) possible to address if it does not render the optimization infeasible.These constraints will need to be satisfied for each time interval t within the optimization horizon in case a time-horizon optimization is concerned based on probabilistic predictions.Besides the constraints on the system states, further time-dependent appliance specific constraints for ∆P t m,q can be included.For example, this has been done in Reference [11].The general optimization function above can be changed to specific optimization objectives depending on the goal and purposes of the DSM application.We will illustrate this by using two examples where each example can be solved using quadratic programming.
Possible objective 1: For the optimization objective of local supply and demand matching over the time-horizon from t = vF to t = vT within each of the households connected to the distribution network while spreading out the remaining supply and demand over time (profile flattening), one can minimize the sum of squares of the base load power vector P t calculated with the flexibility vector ∆P t of each of the prosumers.This is shown in the formula below.min If one alternatively wishes to balance the supply and demand of the distribution network as a whole, one can minimize the square of the sums of the base load power injection P t calculated with the flexibility ∆P t (i.e., not the sum of squares but the square of sums).
Possible objective 2: In a similar fashion, with the objective of fair power sharing, the flexibility provided by different prosumers is aimed to be mostly equal among those who can reasonably contribute to the OLV at hand, which prevents some prosumers from being asked to provide flexibility more often than others depending on their physical point of connection.Since this effect can obviously not completely rule out the physical aspects of the network, some users will be able to make a larger contribution than others (depending on the phase and location of their connection).Therefore, the least sum of squares of ∆P t n,p is taken as a balanced approximation to realize fair power sharing by automatically selecting those prosumers that can reasonably contribute to resolving the OLV at hand.min

Simulation Results
This section aims to demonstrate the effectiveness of the probabilistic time-horizon DSM methodology featuring fair power sharing, which is seen in Equation ( 22) as the optimization objective.The DSM is applied to a three-phase implementation of the IEEE European LV test network as displayed in Figure 4.This network has a radial architecture consisting of 117 nodes and 55 prosumer households with single phase connections.

Simulation Results
This section aims to demonstrate the effectiveness of the probabilistic time-horizon DSM methodology featuring fair power sharing, which is seen in Equation ( 22) as the optimization objective.The DSM is applied to a three-phase implementation of the IEEE European LV test network as displayed in Figure 4.This network has a radial architecture consisting of 117 nodes and 55 prosumer households with single phase connections.

Overall Simulations
Within the simulation setup, the slack node voltage is assumed to remain constant at 230 V and the DSO is assumed to have set voltage limits of 0, 9, and 1 p.u. (i.e., +/−23 V) for the whole feeder.
The acceptable probability limit  for operational limit violations of the voltage magnitudes is set to 0.1 for each time interval (i.e., 90% certainty of having no operation limit violation).Additional constraints on operational limits can be considered such as constraints for the branch current magnitude or the unbalance factor (VUF) for the voltage between the three phases.However, in this work, we focus on violations of the voltage magnitudes in either of the three phases.The assumption is that flexibility is available at each household with an upper bound of 2 kW where future work will be done to make this flexibility more realistic, which is discussed in Reference [37].It should be noted that requests for that much flexibility will only be met in the case of very high loadings of the network such as in cases of multiple charging electrical vehicles.In those cases, 2 kW is less than a quarter of the total household loading and, therefore, considered realistic.For now, the goal is to model flexibility as realistic as possible but to show the effectiveness of the proposed probabilistic sensitivity method.The overall simulation consists of the following steps, which corresponds with the flow

Overall Simulations
Within the simulation setup, the slack node voltage is assumed to remain constant at 230 V and the DSO is assumed to have set voltage limits of 0, 9, and 1 p.u. (i.e., +/−23 V) for the whole feeder.
The acceptable probability limit f for operational limit violations of the voltage magnitudes is set to 0.1 for each time interval (i.e., 90% certainty of having no operation limit violation).Additional constraints on operational limits can be considered such as constraints for the branch current magnitude or the unbalance factor (VUF) for the voltage between the three phases.However, in this work, we focus on violations of the voltage magnitudes in either of the three phases.The assumption is that flexibility is available at each household with an upper bound of 2 kW where future work will be done to make this flexibility more realistic, which is discussed in Reference [37].It should be noted that requests for that much flexibility will only be met in the case of very high loadings of the network such as in cases of multiple charging electrical vehicles.In those cases, 2 kW is less than a quarter of the total household loading and, therefore, considered realistic.For now, the goal is to model flexibility as realistic as possible but to show the effectiveness of the proposed probabilistic sensitivity method.The overall simulation consists of the following steps, which corresponds with the flow diagram in Figure 1.
Step 1: As a first step, the DSO performs the probabilistic load prediction to form the input distributions for the grid supportive demand side management.In the simulations performed in this work, the PDFs are normally distributed [38] where the mean and variance is different for each household and time interval and fitted from historical data obtained from the Pecan Street project [39].The mean values of the household consumption PDFs range from −5.3 to 6.8 kW depending on the household and time of the day while the variance goes up to 5 × 10 6 W 2 with an average of 8 × 10 5 W 2 .Besides the ordinary base load, a significant amount of PV generation is present in the load profiles, which is a high peak load in the evening due to electrification of cooking and heating installations.
Step 2a: For benchmarking purposes, the DSO also performs a Monte Carlo based PPF network for the DA 24-h period by using three phase unbalanced power flow calculations [26].The Monte Carlo simulations for each time interval take 10,000 samples from the PDF of the power consumption of the connected households.From the resulting PDFs of the system states (being nodal voltages and branch currents), the DSO will determine whether the probability of operational limit violations occurring in the network is acceptable and, if not, determines the sensitivity operation point for calculation of the Jacobian.
Step 2b: In the machine learning based approach, the DSO will extract the distribution shape parameters (e.g., mean and standard deviation) and perform the NN based analysis of the IEEE European LV test network for the active power flexibility described in Section 4 based on the PDFs for each of the households.From the resulting outputs of the NN, the DSO will determine for each time interval whether the probability of any operation limit violation occurring in the network is acceptable or not.
Step 3: If the probability for operation limit violations is too high based on the results of step 2a and 2b, the DSO will determine the active power sensitivity of the OLVs at hand.
Step 4: The results of the analysis for OLVs and their active power sensitivities are sent to the DSM application.
Step 5: The DSM application optimizes the flexibility provided by the prosumers in accordance with its own optimization objective.As stated, the optimization objective featured in this work is fair power sharing among the different prosumers that can reasonably contribute in resolving the OLV at hand, which is specified in Equation (22).This way, the optimization will resolve the OLV expected in the network while dividing the burden for doing so over the different participating prosumers.For each appliance in the optimization, additional constraints can be set such as described in Reference [11].
Step 6: As a final step, the effectiveness of the allocated flexibility is assessed.The allocated flexible power comes on top of the original base load, which was represented by the PDFs of the power consumption of the households.Since this base load will still have the associated uncertainty after allocation of the flexibility, the resulting ∆P t of the optimization is added to the original power values of the PDFs of the household consumptions.After this, there is a final verification in the simulation for both the PPF approach as well as the NN approach.It should be noted that the final verification carried out in this work is merely to verify the performance of the proposed approach and will not be carried out during a practical application.During the final verification, for both the PPF approach as well as the NN approach, a Monte Carlo simulation is carried out with the only difference that the input samples from the PDFs are now shifted over ∆P t where ∆P t is either the result of the PPF approach or the NN approach.The overall results of the simulation are discussed in the next section by starting off with the benchmarking results.

Benchmarking Results Using Probabilistic Power Flow
For each of the Monte Carlo simulations within the benchmarking PPF approach, the minimum respectively maximum values out of the vectors [V t,l a ; V t,l b ; V t,l c ] and [V t,u a ; V t,u b ; V t,u c ] are determined for each time interval t and their probabilities and are displayed in Figures 5 and 6, respectively.In other words, Figure 5 displays the probability of the lowest voltage, according to any node or phase within the network.In addition, shown is the probability of having a lower or equal voltage of exactly f .That means that there is a 90% chance that there will be no lower voltage than the displayed value anywhere in the network.If the displayed value is lower than the voltage limit, there is a higher than 10% chance for under voltages.
Energies 2018, 11, x FOR PEER REVIEW 13 of 17 higher than the voltage limit).The dashed blue lines represent the voltage magnitudes over time as they would be when no DSM is applied while the solid blue lines represent the voltage magnitudes after application of the fair power sharing DSM optimization.It should be noted that the displayed values can relate to any node and phase in the network.This can concern a different node and phase at each time interval and even a different node and phase before and after DSM.It cannot be seen from the graphs which node is concerned since this is changing over time depending on the loading of the network and at which node and phase the most extreme voltages (i.e., lowest or highest voltage) occur.From the figures, we can derive that the network is expected not to be capable of facilitating the large PV infeed on a sunny day as installed for this configuration nor the energy intensive appliances that consume power in the evening since there is a higher than 10% probability of having over voltages during the middle of the day and under voltages during peak hours in the evening.During the remainder of the day, there is a less than 10% probability of having operation limit violations, which is where the displayed voltages are below the limit.After application of the DSM optimization, all OLVs have been (nearly) resolved by reducing the probability of operational limit violations for the voltage magnitude below 10% at nearly all times.Some deviations and small OLVs do remain, which can mainly be attributed to the linearization process and the probabilistic uncertainty in the PPF.Energies 2018, 11, x FOR PEER REVIEW 13 of 17 higher than the voltage limit).The dashed blue lines represent the voltage magnitudes over time as they would be when no DSM is applied while the solid blue lines represent the voltage magnitudes after application of the fair power sharing DSM optimization.It should be noted that the displayed values can relate to any node and phase in the network.This can concern a different node and phase at each time interval and even a different node and phase before and after DSM.It cannot be seen from the graphs which node is concerned since this is changing over time depending on the loading of the network and at which node and phase the most extreme voltages (i.e., lowest or highest voltage) occur.From the figures, we can derive that the network is expected not to be capable of facilitating the large PV infeed on a sunny day as installed for this configuration nor the energy intensive appliances that consume power in the evening since there is a higher than 10% probability of having over voltages during the middle of the day and under voltages during peak hours in the evening.During the remainder of the day, there is a less than 10% probability of having operation limit violations, which is where the displayed voltages are below the limit.After application of the DSM optimization, all OLVs have been (nearly) resolved by reducing the probability of operational limit violations for the voltage magnitude below 10% at nearly all times.Some deviations and small OLVs do remain, which can mainly be attributed to the linearization process and the probabilistic uncertainty in the PPF.Similarly, Figure 6 displays the probability of the highest voltage of any node or phase including the probability of having a higher or equal voltage of exactly f (i.e., there is a 90% chance that there will be no higher voltage and a higher than 10% chance for overvoltages if the displayed value is higher than the voltage limit).The dashed blue lines represent the voltage magnitudes over time as they would be when no DSM is applied while the solid blue lines represent the voltage magnitudes after application of the fair power sharing DSM optimization.It should be noted that the displayed values can relate to any node and phase in the network.This can concern a different node and phase at each time interval and even a different node and phase before and after DSM.It cannot be seen from the graphs which node is concerned since this is changing over time depending on the loading of the network and at which node and phase the most extreme voltages (i.e., lowest or highest voltage) occur.
From the figures, we can derive that the network is expected not to be capable of facilitating the large PV infeed on a sunny day as installed for this configuration nor the energy intensive appliances that consume power in the evening since there is a higher than 10% probability of having over voltages during the middle of the day and under voltages during peak hours in the evening.During the remainder of the day, there is a less than 10% probability of having operation limit violations, which is where the displayed voltages are below the limit.After application of the DSM optimization, all OLVs have been (nearly) resolved by reducing the probability of operational limit violations for the voltage magnitude below 10% at nearly all times.Some deviations and small OLVs do remain, which can mainly be attributed to the linearization process and the probabilistic uncertainty in the PPF.
Deviations from the linearization process occur especially when the network becomes highly unbalanced, which is the case during the PV infeed.Nevertheless, the DSM application is slightly conservative most of the time, which overcompensates more when there is a more severe operational limit violation.This can be considered as favourable, as with this, the system operates on the safe side.In the rare event that a small operation limit violation will remain after the time-horizon DSM, additional corrective real-time DSM can be triggered based on state estimation.

Results Using the Neural Network
In the previous subsection, from the benchmarking results, the PPF based approach has been shown to be effective for reducing the probability of OLV to acceptable levels on a DA basis.In this case, the performance of the NN based approach is discussed where the results are presented in Figure 7 for under voltages and in Figure 8 for over voltages.This time, only the voltage corresponding with a probability exactly f is shown.Similar to Figures 5 and 6, as soon as this voltage crosses the indicated voltage limits, the probability of having a OLV higher than f and flexibility needs to be procured.In the figures, the dashed blue lines represent the voltages before DSM while the solid blue lines represent the voltages after DSM using the PPF approach.Lastly, the red lines are for the voltages after DSM using the NN based approach.The insets give a detailed comparison between both approaches during the hours at which operation limit violations take place.Deviations from the linearization process occur especially when the network becomes highly unbalanced, which is the case during the PV infeed.Nevertheless, the DSM application is slightly conservative most of the time, which overcompensates more when there is a more severe operational limit violation.This can be considered as favourable, as with this, the system operates on the safe side.In the rare event that a small operation limit violation will remain after the time-horizon DSM, additional corrective real-time DSM can be triggered based on state estimation.

Results Using the Neural Network
In the previous subsection, from the benchmarking results, the PPF based approach has been shown to be effective for reducing the probability of OLV to acceptable levels on a DA basis.In this case, the performance of the NN based approach is discussed where the results are presented in Figure 7 for under voltages and in Figure 8 for over voltages.This time, only the voltage corresponding with a probability exactly  is shown.Similar to Figures 5 and 6, as soon as this voltage crosses the indicated voltage limits, the probability of having a OLV higher than  and flexibility needs to be procured.In the figures, the dashed blue lines represent the voltages before DSM while the solid blue lines represent the voltages after DSM using the PPF approach.Lastly, the red lines are for the voltages after DSM using the NN based approach.The insets give a detailed comparison between both approaches during the hours at which operation limit violations take place.From the figures, it can be seen that the NN-based approach resembles the performance of the PPF approach in a good way and, in this simulation, it is especially accurate for the over voltages.During the evening hours where under voltages occur, in the second half of the under voltage period, a situation occurs in which the deviation from the PPF approach is higher (although more accurately to the operation limit).In this case, the occurring loading configuration was insufficiently present in the training data for the NN, which results in the lower accuracy.

Execution and Training Time
As has been discussed, the PPF approach is computationally intensive despite the many works in literature improving its efficiency.Carrying out DA PPF calculations for its complete distribution network, this might require significant computational power.Notwithstanding the large computation time for the PPF approach, one may argue that, since the time-horizon preventive DSM proposed in this work runs on a DA/ID basis and not in real-time, computation time is not a top priority.However, as predictions on the load probability density functions tend to be more accurate closer to the time of delivery/consumption, deferring the procurement of flexibility by the DSO as much as possible is important.The DSM application itself should be allowed time to perform its optimization.For the numerical simulations performed in this work, the PPF and DSM optimization for the 1440 time intervals takes over two hours on an Intel Core i7-4790 CPU excluding the verification step in step 4. The computation time is strongly dependent on the number of time intervals for which an OLV is expected.For the NN-based approach, the computation time reduces strongly to only under 15 min because of the fast evaluation of the NN.Besides the execution time, the NN also requires training time.By training the NN on a Nvidia GTX 1080 Ti graphics card, within one hour, training the mean squared error of the voltage magnitudes drops below 0.12 V 2 by using the mentioned gradient descent algorithm.It should be mentioned, however, that the choice for the number of layers and neurons for a particular network configuration might require repeated retraining of networks.Lastly, preparing the training data set might require way more time than the training itself.Still, each of these tasks are in principle a one-time exercise that will save considerably on the computation time during operation of the presented grid-supportive DSM.

Recommendations and Conclusions
This paper contributes with a probabilistic approach for grid supportive demand side management by presenting a Monte Carlo as well as a Neural Network based approach to reduce the probability of geographical dependent operation limit violations to acceptable levels.A practical case study is presented using the IEEE EU LV test feeder for which numerical simulations show that the proposed approach forms an effective method for DSOs to invoke time-horizon flexibility to address (too high probabilities of) expected operational limit violations in the network.In this case, the NN-based approach offers a significant benefit over the PPF-based approach in terms of computational complexity.Nevertheless, from these findings, this research can be extended in several directions.First, the proposed approach will need to be extended with models of actual flexible appliances to allow for more realistic modelling of the available appliance flexibility.Furthermore, the currently used constraints to prevent violation of the nodal voltage magnitudes can be extended to various other power quality related limits such as for the voltage unbalance factor between the three phases.Lastly, a more in-depth quantization of the severity of operation limit violations that can be expected in the future is required together with an analysis of whether the expected available flexibility will be enough to reduce the probability of these operational limit violations sufficiently.The proposed approach needs to be extended to deal with a shortage of flexibility by rendering the optimization infeasible.As a future solution, the optimization problem could be reformulated on the fly such that it will minimize the probability of operational limit violations.

Figure 1 .
Figure 1.Comparison of the proposed grid-supportive DSM using PPF and machine learning.

Figure 1 .
Figure 1.Comparison of the proposed grid-supportive DSM using PPF and machine learning.

Figure 2 .
Figure 2. Illustration reducing probability of under voltages to acceptable levels.

Figure 2 .
Figure 2. Illustration reducing probability of under voltages to acceptable levels.

Figure 3 .
Figure 3. Schematic overview of the separated NN for the nodal voltage magnitudes and angles.

Figure 3 .
Figure 3. Schematic overview of the separated NN for the nodal voltage magnitudes and angles.

Figure 6 .
Figure 6.90% probability maximum voltage magnitudes over time.

Figure 6 .
Figure 6.90% probability maximum voltage magnitudes over time.

Figure 6 .
Figure 6.90% probability maximum voltage magnitudes over time.