Distribution Locational Marginal Pricing (DLMP) for Unbalanced Three-Phase Networks

This paper applies the principles of distribution locational marginal pricing (DLMP) to unbalanced three-phase distribution networks. We first propose a linear model for AC optimal power flow derived through a series of approximation and reformulation techniques. Then a scenario-generation algorithm is proposed to properly model the uncertain parameters in the linear model. Through a proposed No U-Turn sampler (NUTS) based algorithm, probability density functions (PDFs) of DLMPs are calculated. These PDFs provide statistical information about the locational and temporal price risks. By means of applications to two IEEE unbalanced test networks, the numerical results show promising performance for the proposed linear model and the NUTS-based algorithm in creating PDFs of DLMPs. DLMP price densities will be increasingly useful as distribution system operators seek flexible, low risk solutions from embedded generators and aggregators of distributed energy resources.


Constants
Relative cost of segment m in cost curve of DER g; C Marginal operational cost of charging/discharging PV unit v with BSS ($/MWh); C  O NE of the most remarkable changes in the operations of the power system supply chain has been in the distribution networks. For many years they were operated in a passive "invest and forget" manner, being the necessary infrastructure carriers between the high voltage grids and the end-users who simply consumed. But with end-users in the low voltage distribution networks now having embedded generation facilities, both conventional and renewable, with electric vehicles and batteries actively charging and discharging, as well as smart energy management systems, the power flows within the distribution networks exhibit higher volatility and greater quality fluctuations. Furthermore, even the use of the existing infrastructure is changing with various, possibly surprising, innovations. Thus, for example in London, the existing street lamps have been converted into electric vehicle (EV) charging pods [1], with all the awkward implications for the pre-existing network which that entails. As a consequence of all of these end-user changes, distribution power flows are becoming much more challenging to manage.
One traditional, but quite unsatisfactory, solution would simply be to over-invest. However, the direction of the change is now much harder to anticipate and thus any selective network strengthening through fixed investments is risky. There is clearly a need to be more efficient in managing the distribution assets, but optimization models to support this have to be very detailed in terms of their spatial-temporal operations, as well as stochastic in order to represent the random effects of wind, solar, and EV in particular. Furthermore, the value of better optimization models is not just in more precise monitoring of the potential risks in power quality. Distribution system operators (DSOs) are increasingly avoiding the costs and risks of over-investment by procuring flexibility services from the distributed resources [2]. This requires them to accept offers from the owners of embedded facilities for their services in relieving congestion at specific times and locations. Computing a fair value for these services, in terms of economic theory, requires an application of distribution locational marginal pricing (DLMP), analogous to the nodal locational marginal price (LMP) which has become a standard approach for congestion in the wholesale market at the transmission grid level.
The complexity of the analysis is enhanced not only by the large-scale nature of the problem, when the required granularity for all voltages and flows is specified, but also because distribution networks are potentially highly unbalanced due to single/double-phase lines and loads. Thus, DLMP models need to be formulated for three-phase unbalanced distribution networks (TUDNs) to be realistic for the existing distribution networks. Therefore, the computational challenges of DLMP are substantial. Currently there appears to be no practical approach to calculate DLMPs, with all the complexities in real distribution networks. In the absence of even a good approximation to DLMP, the pragmatic choice is often long run marginal cost analysis (minimum increase in total cost when more electricity is consumed by one customer) [3]. This analytical inadequacy creates inefficient signals for consumer behavior, network investment, and the financial performance of distribution utilities [4], [5]. Nevertheless, even without a DLMP solution, utilities recognize a need to adapt their pricing schemes and are motivated, in particular, by the charging and discharging behavior of EVs. For example, 29 distributors in New Zealand planned new pricing solutions (e.g. based on long run marginal costs) in 2020-21 [6], [7]. The ideal DLMP approach would however provide locational and time-of-day prices to keep overall costs down, maintain power quality and avoid the distributional effects of some consumers effectively subsidizing the behavior of others. In this paper we provide such an approach based upon DLMPs, demonstrated on two case-studies to be scalable for real applications.
In order to represent the stochastic effects, e.g. from wind and solar outputs, on the DLMPs, probability density functions (PDFs) of DLMPs are estimated. The probability density function f X (X) of a DLMP, specified on the continuous random variable X, can be used to estimate the extreme local price risks. As regulated companies, DSOs will be cost-conscious and will be averse to the risk of procuring flexibility services at high prices. The conventional risk management controls in many companies use the so called "value-at-risk" limits, computed from the PDFs. Thus a 95% value-at-risk limit a would be the value at which there is only a 5 percent probability of the DLMP being higher. If b is the market price cap, it would be expressed as Therefore, these PDFs provide information regarding the risks of high or low prices in different parts of distribution network, being useful, we would argue, not just to the risk-averse procurement of services by the DSO, but also to the asset owners of flexibility services in terms of their operational revenue risks. Some researchers such as [8] have gone beyond point estimates of the mean nodal prices to include variance as a risk measure, but this is clearly inadequate to properly estimate the tail probability risks in non-Normal distributions. Hence, we propose to use the full PDFs. Overall, such probabilistic analysis through PDFs provides a more complete assessment of TUDN risks for both operational and investment decisions.

B. Background Research
Power flow equations are intrinsically nonlinear (due to the convolution of voltages) and so determining DLMP becomes a  nonconvex and nonlinear problem. Whilst various approaches can handle this nonlinearity, we have focused on developing a linear program (LP) approximation. From existing research, second order cone programming (SOCP) [9], semi-definite programming (SDP) [10], and branch flow formulations [11] are the most common alternatives. The algorithmic tightness and solution times are crucial considerations in the different solution algorithms [12]. Nevertheless, these approaches may be employed for three-phase power flow models [13]. Regarding DLMPs, these are mostly calculated for balanced distribution networks such as [14] which uses the branch flow formulation. But distribution networks are naturally highly unbalanced, since the majority of consumers are single phase demands. Although Hanif et al. formulated a multiphase DLMP model based on SOCP [15], it requires heuristics to deal with the tightness and solution times at large-scale. Another class of methods to determine LMPs is predictive, based upon historical data. This is an active research theme in which the prices are predicted either by 1) applying machine learning approaches such as neural networks [16] or 2) simulating the network and market behavior [17], [18].
Traditionally researchers develop point estimates for the LMPs [19], [20] without addressing the uncertainties in the prices. However, with the rise in intermittent end-user generation, these uncertainties are becoming more important operational considerations. Thus, machine learning and maximum likelihood methods have been used to forecast the variances as well as means of the prices [21]. We have summarised several relevant papers in Table I, in which it is evident that most of the research focus is on single phase models assuming that the distribution network is balanced. The only available three-phase models for distribution networks can be broadly categorized as: (1) A group of models proposing an SDP for modeling three-phase distribution networks. The SDP models are computationally complex and it is often hard to find a solution with zero optimality gap. (2) The second group of models rely on DC optimal power flow (DCOPF) approximation. While this approximation might be acceptable for transmission network operation, it is not suitable for distribution networks due to their high R/X ratio. The R/X ratio in a distribution network is between 0.5 and 2 while it is less than 0.1 in a transmission network.
(3) The third group of models solve the ACOPF formulation for distribution networks using variants of the Newton-Raphson (NR) method. However, the NR method has a slow convergence rate and its convergence to a solution is not guaranteed.

C. Contributions of This Paper
This paper contributes to the existing literature in the following ways: 1) We propose a linear mathematical model for the operation of three-phase unbalanced distribution networks (TUDNs). To the best of our knowledge, the proposed model in this paper for calculating DLMPs is state-of-theart. It is applicable to three-phase unbalanced distribution networks and it does not suffer from the convergence issues, relaxation gap, and high R/X ratios, as do the various alternative approaches listed in Table I  studies and numerical results and Section IV concludes the paper.

II. THE LINEAR PROGRAMMING MODEL AND THE PROPOSED ALGORITHM
The network optimization problem for a TUDN is formulated as an LP problem in (1).
cφts , v bts , u bts }. The objective function (1a) represents the expected operational cost of the TUDN over different scenarios s ∈ S with probability PR s . The first term is the operational cost of the distributed energy resources (DER) units with fuel cost C (G) g for generator g. The fuel consumption of DER units is formulated as a second order quadratic function of the generated power. Then, this function is linearized with the piecewise-linearization technique as shown in Fig. 1.
mg is assigned to the cost-curve segment m of the DER unit g. The second term models cost of the received energy from the distribution feeder f . The next term is used to penalize the curtailed demand d with a value of lost load C (D) dt . The last terms are marginal costs of charging/discharging different storage systems multiplied by charged/discharged active power. The active power generation of DER unit g at each segment m is constrained in (2b). The total active power generation is constrained in (2a). Since we need the power flow equation for each phase φ ∈ P , the variable p (P ) gφts is defined to model the DER generation of each phase. The total active power generated by DER unit g is the sum of power on all phases as well as all segments as shown in (2b). Similarly, the generated active and reactive power by DER unit g (both emergency and DER units) and distribution feeder f are constrained in (2c).
gφts ≤QP φg UG gts (2c) The emergency generation units have black start capability and could deliver energy even when their buses are not energized. In contrast, other DER units are able to generate electricity only when their corresponding buses are energized. The feeders apparent power (p fφ which is linearized employing the circular constraint linearization method shown in Fig. (1). Accuracy of this approximation technique is discussed in [33]. Eight linear constraints approximate the feeder power capacity constraint in (3).
Similarly, the active and reactive power flowing through the feeder (b, k) ∈ L is limited by (4a)-(4d). Two square constraints, first square constraints are (4a) to (4b) and second one are (4c) to (4d), are employed to substitute (p The binary parameter UL φbkt represents the feeder status which is zero if the corresponding phase is not available (e.g. due to a fault).
The active/reactive demands are controllable and curtailable (UD dt is availability of demand). The binary parameter UD dt models the on/off status of the demand unit. Ohms law for each distribution feeder can be written as v bts = v kts -Z bk i bkt [13]. When both sides are multiplied by their conjugate transpose (The sign () is used for element by element complex conjugate for matrices), we have The last term is nonlinear which will be linearized in the following. When phase voltages are nearly balanced and feeder losses are relatively small Z bk i bkt S bk , the linear form of the unbalanced three-phase power-flow equation can be written as u bts -u kts + Z bk S bk + Z bk S bk = 0 where the terms Z bk i bkt are removed. The vector variable of square-voltage magnitudes u bts = [u bφts ] φ∈P is used to linearize the quadratic term. Similarly, the vector of currents is defined as i bkt = [i bφkt ] φ∈P . The linear power flow model for a TUDN is presented in (5). The binary parameter UL φbkt and the big-M are used to model feeders which are not available at time t.
The symbol is used for component-wise multiplication of matrices. The voltage magnitude is limited in (6) and fixed for slack buses.
The active and reactive power balance is written in (7a) and (7b) where the E b is set of all battery storage systems (BSSs) connected to bus b. The Lagrange dual variables λ bφts /μ bφts of the active/reactive power balance equation in (7a)/(7b) represents DLMP for active/reactive power, i.e. DLMP-P/DLMP-Q.
The reactive power of a shut capacitor c ∈ C is calculated by multiplying square of its voltage magnitude u bφts by its per-phase susceptance B (C) c in (8) which is then constrained by its nominal reactive power. The number of switched steps SC ct of the switchable shunt capacitor c ∈ C is used to calculate reactance of the shunt capacitor bank. Each switchable shunt capacitor c ∈ C has a susceptance of fixed capacitor B (C) c which can be connected to the system all the time, and switchable capacitors which are switched to meet the TUDN requirements. By connecting SC ct steps of switchable capacitors with susceptance B (SC) c , total susceptance of the switchable capacitor will be B The active power output of PV unit v is constrained in (9) to nominal output PV vφ and forecasted output FV vφt based on forecasted solar irradiance. Also, the reactive power output is constrained in (9) by its nominal output QV vφ .
The prosumers with PV units and BSS v ∈ V (E) can store the absorbed energy in their BSS as modeled in (10). Therefore, these units could be used to generate electricity as a grid-forming unit with black start capability during blackout. vφts is limited in (10b). The active and reactive powers of wind unit w ∈ W are limited in (11).
The prosumers with wind units and BSS (w ∈ W (E) ) can store extra generation as modeled in (12). Therefore, these units, also, could be used to generate electricity as a grid-forming unit with black start capability following a blackout.
The set of constraints (13) models a general BSS.
The voltage regulator or tap changing transformer r ∈ R is modeled as three single-phase voltage regulators connected between buses b and k with maximum/minimum ratio vector as R r /R r . Each single-phase voltage regulator is modeled with an ideal transformer in series with a leakage impedance. Therefore, in a three-phase voltage regulator R r R r ≤ u bts /u kts ≤ R r R r . Which is enforced by (14) using the binary on-off parameter of the regulator UR rts and the availability vector A

A. The Proposed NUTS-Based Algorithm to Estimate PDF of DLMPs
The Hamiltonian Monte Carlo (HMC) algorithm is a Markov Chain Monte Carlo (MCMC) method which can be used to propose a sequence of samples that follow a target distribution for which a direct sampling is difficult [34]. The HMC has been used in reliability analysis in [35] and for estimating genetic parameters and breeding values in [36]. As compared to many MCMC methods, the HMC algorithm avoids (a) the randomwalk behavior and (b) sensitivity to correlated random variables. Accordingly the HMC algorithm can effectively explore the target probability space [36]. When estimating a random variable x with probability density function f (x) using the HMC algorithm, we define an auxiliary momentum variable ρ and consider the joint probability f (x, ρ) = f (ρ|x)f (x). In practice, ρ follows a normal distribution f (ρ) ∼ N (0, M) (where M is the covariance matrix). Now, we define the Hamiltonian as H(x, ρ) = − log f (x, ρ) which can be re-written as: . We can solve the Hamiltonian dynamics (equations) below to create a new sample (x, ρ) from an existing one (which is a system of differential equations) [37].
Where t is a fictitious time. The standard approach in the literature is to solve the Hamiltonian dynamics (15) in a discrete time setting using the leapfrog algorithm (there is no analytical solution to solve this dynamic system) [34]. The HMC algorithm at each iteration t solves the Hamiltonian dynamics using discrete approximation to find the new sample at next iteration t + Δt. Fig. 2 shows an example probability space where the x-axis represents the random variable x, the y-axis is the auxiliary momentum variable ρ, and the z-axis (shown with coloring) is the joint probability function f (x, ρ). In Fig. 2, performance of the leapfrog (random-walk) algorithm in two different initial samples (a) and (b) are compared with the HMC algorithm in (c) and (d). To estimate the joint probability function, the HMC starts with an initial sample and solves the Hamiltonian dynamics (15) to find the next sample. As we can see the marked HMC paths follow our target distribution and avoid the random-walk behavior.
However, performance of the HMC algorithm is highly sensitive to proper tuning of the step-size (Δt) and number of steps (NP ) [38] in solving the Hamiltonian dynamics (15). Recently, the NUTS algorithm has been introduced to improve the HMC algorithm [39]. There is no need to tune these parameters (Δt and NP ) in the NUTS algorithm. In this paper, we have used the NUTS algorithm as an extension of the HMC algorithm to estimate PDF of DLMPs. General steps of the proposed NUTS-based algorithm are shown in Fig. 3. It uses the provided data (historical data of WT, PV, demands, and electricity price) as  input and provides PDF of DLMPs as output employing the proposed NUTS-based algorithm. The algorithm starts with gathering the input data. Then, a cluster-based scenario generation is used to generate required scenarios to solve the stochastic programming problem (1). Optimal DLMP-P/Q prices are obtained by solving (1) and finding the Lagrange multipliers associated to active-power and reactive-power balance constraints. At this stage, the initial samples for DLMP-P/Q prices and momentum variables ρ are generated. Then, the NUTS algorithm is used to solve the Hamiltonian dynamics (equations). This algorithm keeps moving to a new point in the probability space to estimate PDF of DLMPs.
No U-turn sampler: In the following, Steps 4 to 6 of the proposed NUTS-based algorithm flowchart are explained in detail. In the proposed NUTS-based sampler (Fig. 4), a multimodal distribution is used as the target distribution for DLMPs (P). This choice is motivated by our observations on the calculated DLMPs in different TUDNs. The total number of iterations (NP ), number of iterations to adapt (NP a ), initial position (x 0 ), and desired mean acceptable probability (δ) are inputs of the algorithm where NP a iterations are used to adapt the step size ξ. A heuristic approach is employed in Steps (iii) to (vi) to find an initial value for step size ξ where the step size is regularly divided by two or doubled until the acceptance probability is reached. Hence, NP iterations are done in Steps (vii) to (xx). At each iteration, a random direction is taken (Step (x)) and this process builds a binary tree via repeated doubling (Steps (xi) to (xiv)) with the proposed tree algorithm in Fig. 5. Doubling stops if the probability of the previous state is very low in Step (xv) or when states overlap (making a U-turn) in Step (xvi). The step size ξ is only tuned in adaptive iterations in Steps (xvii) to (xviii). After NP a iterations, tuning stops and the algorithm keeps using the tuned step size ξ NP a . The proposed NUTS-based algorithm enables us to estimate the PDFs of DLMPs accurately with less computational burden than if all scenarios were considered. The performance of the proposed NUTS-based algorithm is analyzed later in Section III. The proposed NUTS-based algorithm in Fig. 4 uses a tree algorithm in Fig. 5 to build a binary tree via repeated doubling [39]. Each state of the proposed NUTS-based algorithm is described by a position x = [λ bφts , μ bφts ] and a momentum ρ. Inputs are current state (x, ρ), initial state (x 0 , ρ 0 ), resample value (y), direction (z ∈ {−1, 1}), step size (ξ), and tree height (j) which are given to the tree algorithm in Step (ii). The choice of direction z is done with a discrete uniform distribution. For a given real number x, the step function is used as This algorithm uses a leapfrog integrator to discover the next steps by going forward or backward for 1, 2, 4, etc. steps (repeated doubling) via a recursive algorithm. If the tree height is zero in Step (iii), one step is taken in the given direction z in Steps (iv) and (v). If the height jis non-zero, it does recursive doubling with height j − 1. This algorithm stops by zeroing the indicator w when most right state and most left state overlap (making a U-turn). Therefore, if w = 0 in Step (vii) then the tree algorithm returns the calculated values in Step (xi). If the indicator w = 1 in Step (vii), there is no U-turn yet and the leapfrog integrator is used again to take one more step in the given direction z in Steps (viii) to (x). At each doubling, the next forward/backward state (x + /x − , ρ + /ρ − ) is found in Step (viii)/(ix). The previous state x is defined in Step (x). Then the tree algorithm returns the results in Step (xi).

III. CASE STUDIES
In this section, the effectiveness of the above model is evaluated in two TUDNs (modified IEEE 34-and 123-bus unbalanced networks). We compare the calculated PDFs of DLMPs using all scenarios with results of our proposed sampling approach. The feeders and demands in these networks are unbalanced (detailed characteristics are available in [40]). The costs of active and reactive power vary together by time and scenario and they are at most $1/kWh.

A. IEEE 34-Bus Three-Phase Unbalanced Distribution Network
The single-line diagram of the IEEE 34-bus TUDN is shown in Fig. 6. This network is supplied by two DERs, two WTs, two PV units, two BSS units, and the main grid through the point of common coupling (PCC). Their features are summarised in Table II. This network has two three-phase shunt capacitors, two voltage regulators, one voltage transformer, and one of each WT and PV units has a BSS unit. The voltage profile of the 34-bus TUDN is shown in Fig. 7 in p.u. for phases A to C. The PCC at bus 5 is used as the slack bus and its voltage is fixed to 1 p.u. Voltages in the other buses change due to the directions of active/reactive power and losses. This graph shows mean voltage and corresponding 95% confidence interval for each phase. They are within the predefined minimum/maximum of 0.95/1.05 per unit. The voltages drop as low as 0.9551 p.u., due to the relatively   Fig. 7, again with 95% confidence intervals around the means. Generally the flow of feeders corresponds to the changes in demand. The key factors that change over time are shown in Fig. 8. These values are scaled by their maxima to fit in one graph. Each of these daily factor profiles are selected from different days to be representative means for zone SE1 of Sweden. Thus, demand is from 2019-03-01, wind power is from 2018-02-26, solar irradiation is Global Horizontal Irradiation (GHI) from 2019-05-19, and generation cost is from 2019-12-18.
The GHI historical data, taken from the second Modern-Era Retrospective analysis for Research and Applications (MERRA-2) which is a NASA atmospheric reanalysis in [41], are used to calculate the outputs of PV units. Generation cost, wind power, and demand data are from Nord Pool [42]. The uncertainties around these means are modelled with four stochastic parameters, i.e. for wind power generation for the WT unit, generation cost, solar irradiation for PV units, and demand.
The calculated DLMPs represent the marginal effects on total operational cost for a unit change in demand at bus b which vary over time as shown in Fig. 9. Next we investigate the effectiveness of a sampling based approach to estimate these calculated DLMPs, and the first stage of our approach to this is to reduce the size of the problem by means of cluster analysis.

B. Clustering Based Scenario Generation
A clustering based scenario generation method is used to generate the scenarios for wind power, solar energy, demands, and electricity price. First the historical data for these uncertain factors in 2019 (2018 for wind) are normalized. Then the data is sorted into k groups using the k-means clustering algorithm. The so-called silhouette values are calculated to determine the best number of clusters. Thus, in our case, five clusters were selected since 5 has the highest silhouette value in Fig. 9. Means of each factor are calculated at each cluster and then random normal variations are generated accordingly. The number of scenarios is adjusted from 1 to 200 and sensitivity of mean operational cost (total operation cost in $/h divided by number of scenarios) are shown in Fig. 10. With increasing number of scenarios, the total operational cost is gradually stabilizes and we found that there was no need to enhance number of scenarios beyond 20 in this network if we only look into mean total operation cost (but in calculating PDFs we have used 100 scenarios to be sure the results are not sensitive to the number of scenarios).
Next, this network is solved 985 times (with 100 scenarios each time) to see changes in DLMPs. The PDF of DLMPs is shown in Fig. 11 for the 34-bus unbalanced network with a total of 985×100 =98,500 scenarios for simulating a single day. Two buses at different times are shown to illustrate the distributions of prices in three-phases. DLMPs have multimodal distributions which encourages us to use this specification later in the proposed NUTS-based algorithm. Bus 14 is connection of a single-phase (phase B) feeder so other plots are not applicable.
Next we show the value of using the proposed NUTS-based algorithm. This works with a set of correlated samples to estimate DLMPs as a multimodal distribution (defined as a mixture of three normal distributions). Firstly, the convergence of the means and SDs of DLMPs with increasing number of scenarios is analyzed. Fig. 12 shows how the mean and SD of the generated samples by the proposed NUTS-based algorithm converged to the mean and SD of the calculated DLMPs (data). It appears that 300 scenarios (instead of the 98,500 scenarios) by the proposed NUTS-based algorithm is sufficient to match the calculated DLMPs. Secondly, Fig. 12 shows the capability of generating adequate samples by the proposed NUTS-based algorithm to match the calculated DLMP PDFs.
The proposed NUTS-based algorithm is implemented in Python using the PyMC3 package [43] and applied to 300 observed data points (specified by bus and over time) to estimate the multimodal distribution. For DLMP-P and DLMP-Q, the mixing parameters for the three latent normal densities are all non-zero, so the choice of three was appropriate, but four would have been superfluous. As the number of scenarios are decreased from 98,500 to 300, the elapsed time goes down from 14:47 minutes to 3 seconds.

C. IEEE 123-Bus Three-Phase Unbalanced Distribution Network
The IEEE 123-bus TUDN is shown in Fig. 13 (details available in [40]). This network has 118 feeders with 11 types of overhead feeders (one-, two-, and three-phase) and one type of underground three-phase feeder. There are one three-phase shunt capacitor, three single phase capacitors, two transformers, and four voltage regulators.Total active/reactive demands are 1420/775, 915/515, and 1155/635 (10 −3 MW)/(10 −3 MVAr) in phases A, B, and C, respectively.
The voltage profile of this network is shown in Fig. 14. The PCC is used as the slack bus and voltages at the slack bus are fixed to one, whilst the other voltages change due to the feeder flows and losses. This graph shows the voltage means and 95% confidence intervals, which are within the predefined acceptable range 0.95 to 1.05 p.u. Since feeders in this network are relatively smaller (from 30.5 m to 304.8 m with mean 101.2 m and SD 47.9 m), the voltages do not vary as much as in the previous 34-bus TUDN.  This network is supplied by five DERs, four WTs, ten PV units, four independent BSS units, and the main grid through the PCC in bus 610. Detailed characteristics are presented in Table III. Four solar units have BSS as well as three wind generation units. The price of electricity drops during the day from 11:00 to 14:00 due to the output of the solar units. Feeder flows are shown in Fig. 14 for each phase and all of them are within the acceptable range (0.95 to 1.05 p.u.). 500 scenarios are used at each time we solve (1) for the 123-bus TUDN.
Next, problem (1) is solved 60 times (with 500 scenarios each time) to visualize changes in DLMPs. The PDF of DLMPs in the 123-bus TUDN for 60×500 =30,000 scenarios over one day operation time (NT = 24 hours) for buses 13,15, and 59 are shown in Fig. 15 for hours 10:00, 14:00, and 17:00, respectively. As before, the PDFs clearly have multimodal distributions. Buses 15 and 59 are connected to phases C and B respectively. Therefore, the inapplicable phases are shown blank. The prices of electricity do not change with phases. Fig. 16 represents the convergence of the means and SDs of the sampled DLMPs which in turn shows that 500 scenarios are sufficient for obtaining PDFs. The estimated PDFs of DLMPs at bus=15, time=14:00, and phase=C are shown in Fig. 16 for 500 scenarios. Using the proposed NUTS-based algorithm for the 123-bus TUDN reduces the computation time significantly from  617 hours (considering all scenarios) to 10 hours (considering only 500 scenarios using our proposed NUTS-based algorithm). Total number of iterations (NP ), number of iterations to adapt (NP a ), initial position (x 0 ), and desired mean acceptable probability (δ) are set to 4000, 1000, 0, and 0.9, respectively.
In the proposed linear model, the linearization is explained in detail in Section II for the constraints (5), (6), and (7). In this model, the relatively small feeder losses Z bk i bkt S bk are used. Accuracy of the linearization method is validated by calculating | Z bk i bkt S bk | which should be as close as possible to zero for the linearization to be valid. Therefore, value of | Z bk i bkt S bk | is calculated for the IEEE 123-bus TUDN as shown in Table IV. Mean of | Z bk i bkt S bk | for phases A, B, and C are 8.0238 × 10 −6 , 1.0012 × 10 −5 E, and 8.0605 × 10 −6 , respectively. Therefore, | Z bk i bkt S bk | is close to zero and feeder losses are relatively small Z bk i bkt S bk and (5) is an accurate linearization. This evidently shows the accuracy of the proposed model in Section II.

IV. CONCLUSION
In this paper we have calculated DLMPs for three-phase unbalanced distribution networks. As shown in the case studies based upon IEEE unbalanced test networks, DLMPs vary by phase, time, scenario, and location. DLMPs are of increasing interest as distribution network operators seek to manage their networks by making more use of flexibility services from distributed energy resources such as batteries or from aggregators operating several kinds of assets with various response capabilities. These DLMPs can be the basis for the valuation of such services in short term operations and thereby help motivate long-term investment in more flexibility assets and arrangements. Furthermore, with increasing uncertainty in consumer behavior through local generation, storage, and EV activity, the need to model these DLMPs as PDFs is becoming crucial for the network operators. The PDFs would allow a risk analysis and stochastic optimization of the flexibility resources. Risk averse DSOs and owners of these flexibility resources will be concerned about the risk of extremely high or low local prices. The PDFs are necessary to compute these risks. In this paper, we have focused upon an efficient method of generating these PDFs, which may be used operationally in various ways. In particular we demonstrated the effectiveness of a new sampling approach. Following a scenario-based cluster analysis, DLMPs are calculated through a NUTS-based algorithm for active and reactive powers. The numerical results demonstrate the accuracy and computational efficiency of the proposed linear model in two IEEE unbalanced distribution networks. The performance of the proposed approach is remarkable in maintaining accuracy at much lower computational times, and this suggests that the scalability to real systems will be feasible. Overall, we contend that an efficient sampling procedure such as the one developed in this paper is crucial if PDFs are to be used in larger scale practical applications. It is not just the size and complexity of the network that creates this computational requirement, it may also be a matter of timeliness, as such prices may be required at high frequency throughout the day as circumstances on the network change rapidly. In future paper, authors are planned to compare the DLMP results with other tariffs strategies and improve the proposed scenario generation algorithm.

APPENDIX: ILLUSTRATIVE EXAMPLE
An illustrative example is used in this section to explain the proposed NUTS-based algorithm in Fig. 4 for estimating the PDFs of DLMPs. It is 100 sample points from the normal distribution with mean=0 and standard deviation=1 (N (mu actual = 0, sigma actual = 1)). The sample points and their PDF are shown in Fig. 17.
Parameters of the proposed NUTS-based algorithm in Fig. 4 and Fig. 5 are P : mu = N (0, 1), sigma = 1, return = N (mu, sigma), NP = 50, NP a = 5, x 0 = 0, and δ = 0.95.  NP a iterations are used to adapt the step size to ξ = 0.09896. Iterations 1 to 8, from NP iterations, in the proposed NUTSbased algorithm are shown in Fig. 18. Current and proposed means are shown at each iteration.
The PDF of the actual data and one set of the generated samples in the proposed NUTS-based algorithm are shown in Fig. 19 which demonstrates similarity of the PDFs obtained from the actual data and from our NUTS-based algorithm. Mean values are 0.12 and 0.09 in the generated samples and the actual data, respectively. 50 predictive samples are generated employing trace of the model P. The PDF of all these samples is shown in Fig. 20 with black lines and PDF of the actual data is shown with blue  lines. As we can see, the proposed NUTS-based algorithm is able to generate samples with similar probability distribution as the original samples in Fig. 17.
As another numerical study, performance of our proposed NUTS-based algorithm is compared with the leapfrog algorithm. Accordingly, the leapfrog algorithm is also applied to this illustrative example. Iterations in the leapfrog algorithm are shown in Fig. 21. Due to the random-walk feature of leapfrog algorithm, proposed means keep changing at each iteration. Current and proposed means at each iteration are also shown in Fig. 21. The current mean is shown with blue circle and the proposed one is shown with either green circle (if accepted) or red circle (if denied).
After 50 iterations, in both algorithms, PDFs of sample means in the leapfrog and our proposed NUTS-based algorithms are shown in Fig. 22. Mean values are 0.12 and 0.22 in the proposed NUTS-based and the leapfrog algorithms, respectively. Also PDF of the samples in the leapfrog algorithm is different from the one from the NUTS-based algorithm as shown in Fig. 22. This is while our proposed NUTS-based algorithm had similar PDF to the actual data as was shown in Fig. 19.