Shapley value cooperative game theory-based locational marginal price computation for loss and emission reduction

An iterative method based on Shapley Value Cooperative Game Theory is proposed for the calculation of local marginal price (LMP) for each Distributed Generator (DG) bus on a network. The LMP value is determined for each DG on the basis of its contribution to reduce loss and emission reduction, which is assessed using the Shapley Value approach. The proposed approach enables the Distribution Company (DISCO) decision-maker to operate the network optimally in terms of loss and emission. The proposed method is implemented in the Taiwan Power Company distribution network 7 warnings consisting of 84 buses and 11 feeders in the MATLAB environment. The results show that the proposed approach allows DISCO to operate the network on the basis of its priority between the reduction of active power loss and emission in the network


Introduction
Restructuring electrical power utilities, increasing load demand, technical innovation, concern for greenhouse gas emissions, improving voltage and loss reduction are some of the key reasons for increasing the penetration of Distributed Generator (DG) units into the distribution network [1][2][3]. In addition, DG Units can provide customers with freedom and flexibility in planning and developing installations according to load criticality [4]. Following the integration of DGs, distribution networks are transformed from passive to active states, similar to transmission networks. Some of the practices that have been implemented in transmission systems, such as nodal pricing, may also be applied in active distribution networks.
The role of the Distribution Company (DISCO) decision-maker is very important in a deregulated environment, and some of the responsibilities are • Ensuring that the network operates optimally in terms of certain parameters, such as loss and emission. • Ensuring that the decisions taken must be economically viable for both DISCO and DG owners • Maintaining fair competition between owners of private DGs and control of privately owned generators.
The above responsibilities can be fulfilled by providing financial incentives to each DG owner. Nodal pricing is one of the well-organized procedures for financial incentives [5] and various policies are available including the well-know Locational Marginal Price [6,7]. Distribution LMP formulation based on active and reactive power losses is developed in [8], and linearized power flow for distribution (LPF-D), loss factors for distribution (LF-D), and linear optimal power flow for Page 2 of 11 Veeramsetty Protection and Control of Modern Power Systems (2021) 6:33 distribution (LOPF-D) are processed to compute distribution LMP in a distribution system. In [9], three-phase distribution LMP computing mechanism is developed to operate the distribution network efficiently by minimizing the losses, improving the voltage profile and preserving the balance between the three phases. In [10], an extended Kalman Filter based method is proposed to determine the nodal prices at DG buses in the distribution network based on the reduction in active power losses by minimizing the total merchandising surplus (MS). The three methodologies have made important contributions to computing LMPs. However, they have not considered the active distribution system for simultaneous monitoring of loss, emission and reliability. A new methodology is proposed to compute locational marginal price in [11] based on energy price, loss prices caused by nodal active power and reactive power, congestion price and voltage support price using linearized AC power flow. Power-based distribution locational marginal pricing (PDLMP) to determine the active and reactive power tariffs in distribution systems with high penetration of distributed energy resources is proposed in [12], where PDLMP for DG owners is computed based on over-voltages, congestion, renewable energy share and distribution system operation cost. Most publications are on LMP computing in distribution network. The comparison of the research contributions of different authors to the different features considered is shown in Table 1  The method proposed in [17] contains all the mentioned parameters, and the nucleolus theory is used as a solution to the problem of cooperative game theory. The main drawback of the nucleolus theory is that it is not monotonous, and the incentives for each DG owner are compared using wholesale market price ( t ) directly. This form of calculation gives each DG owner more incentives, but reduces the income of the DG owner for some time.
In this paper, the Shapely Value Cooperative Game Theory (SVCGT) based iterative method is proposed. This is an improved version of the method in [17]. Pursuant to the proposed method, the Shapely Value solution concept, which is not monotonic as opposed to nucleolus, is used to find the share of each DG unit in loss and emission reduction. The incentives for each DG are then calculated by distributing the financial savings of DISCO because of loss and emission reductions. This improved method also includes all the features of the comprehensive method proposed in [17]. SVCGT is used for the first time to compute LMP based on active power loss and emission reduction, and the proposed method can benefit the DISCO in the following ways. It: • Can simultaneously run the network with optimal active power loss and emission; • Can sustain fair competition between private DG owners;  • Can monitor private DG owners on the matching between the injected power from the DG units into the network and the computed generation on the basis of priority between loss and emission reduction.
The main contributions of this paper are: • New mathematical modeling is developed to compute the incentives for DG owners; • Incentives for DG owners are calculated on the basis of DISCO's financial savings due to loss and emission reduction: • The Shapely Value method is used to identify the contribution of each DG unit to the reduction of loss and emission.
The problem formulation and methodology are provided in Sect. 2, while the discussion on simulation results and the effectiveness of the proposed method are in Sect. 3, Finally, conclusions are drawn in Sect. 4.

Problem formulation and methodology
The problem addressed in this paper is on the computation of LMP at DG buses in a radial distribution system based on loss and emission reduction, and is approached from the two aspects of : • Allocation of share of each DG unit in loss and emission reduction from the base case using the Shapley method • Calculation of incentives based on the share of each DG unit in loss and emission reduction from the base case.

Load Forecasting
In order to estimate the LMP at each hour of the next day, it is necessary to predict the corresponding load at each period. Two layers of the Artificial Neural Network (ANN) are used to predict the load, with the expected load of L(i, D) considered to be the output, as shown in Fig. 1. As seen, ANN is designed and trained to predict load L(i, D) based on the loads in the previous four hours and at the same hour in the previous two days. However, the forecasted price is considered from [22]. A Back Propagation Algorithm (BPA) is used to train the ANN because of its flexibility and, learning capability, and is highly suitable for problems with no relationship between output and input [23]. The efficiency of the ANN is measured in terms of MAPE and RMSE [24][25][26]. A stochastic gradient descent approach is used to train the model, and the training dataset consists of 6644 samples and the testing dataset consists of 24 samples.

Calculation of loss reduction
To calculate the loss reduction, an efficient load flow solution must be run with a projected load of L(i, D) in two cases as: • Base case: Distribution network without DG • Case2: Distribution network with DG We use the backward and forward load flow solution [27] as it takes full advantage of the ladder structure of the distribution network to achieve high speed, robust convergence and low memory requirements [28,29]. Before running the load flow for Case2, the generation of each DG unit based on the active power price and base case loss reduction are calculated respectively using:

Computation of emission reduction
Emission in the base case is computed based on the load on the substation bus and the emission coefficients corresponding to the substation bus as shown in (3), whereas emission in Case 2 is computed based on the power injected from the substation bus and generators, the emission coefficients for the substation bus and DG units as shown in (4). Emission reduction of case 2 from the base case is shown in (5).
(1) The emission penalty for DISCO in the base case and in Case 2 are calculated using (6) and (7), respectively, while financial saving for DISCO due to emission reduction from base case is computed using (8), i.e.:.

Shapley Value Method
The Shapley Value (SV) method [30] is used to allocate the share of each DG in loss and emission reduction shown in (2) and (8), respectively. The assignment of the share in reduced loss and emission among all DG units is carried out using: where, for reduced loss allocation, k=l, and for reduced emission allocation, k=e.
The allocation of the share of each DG in loss reduction can be explained by sample calculations. Considering a distribution system with three distributed generators and a basic case loss of 220 kW, Table 2 shows the loss reduction due to each DG unit coalition as shown in [18].
The share of each DG unit in loss reduction computed as in (9) is shown below: As seen, shares of the DGs in loss reduction when all units participated in coalition are 33.0kW, 51.9kW and 27.2kW, respectively. This leads to a total loss reduction of 112.1kW as shown in Table 2. A similar share of each DG in emission reduction can also be computed using the SV method.

DISCO's Extra Benefit
DISCO's extra benefit is defined as the difference between DISCO's benefit with and without DG units. DISCO's benefit without DG is defined as the difference between the amount paid to obtain power from the main grid and the emission penalty and the revenue collected from customers. Similarly, DISCO's benefit with DG is the difference between the amount paid to obtain power from the main grid and DG units, the emission penalty and the revenue collected from customers. Thus, DIS-CO's benefit with DG, DISCO's benefit without DG and DISCO's extra benefit can be calculated respectively as: Similarly, DISCO's extra benefit in the presence of DG units due to loss and emission reduction can be minimized by providing incentives to DG owners.
In the case that DISCO provides power to customers from low emission coefficient generators, customers will pay a lower energy price because of a lower emission penalty. Similarly if customers receive power from high emission generators a higher energy price will have to be paid. As both customers and generators must share emission penalty, the emission penalty sharing factor ω e is considered as 0.5 here in this paper, indicating equal sharing between the generators and customers.

SVCGT based iterative Algorithm
The computation of LMP at each DG bus using the Shapley Value method, based on loss and emission reduction is explained in the iterative algorithm. (12)     The proposed SVCGT based iterative method for LMP computation at DG buses is simulated on the Taiwan Power Company (TPC) distribution system consisting of 84 buses and 15 DG units. Complete data of the TPC distribution network has been considered from [32]. Cost coefficients of each type of DG are shown in Table 3., and emission coefficients are shown in Table 4.
With the assumption that all 15 DG units with 1MW capacity are operated with 0.9 leading power factor and are placed in the network as shown in Table 5, the proposed algorithm is simulated in MATLAB [33] on a i5-5200U computer having 2.25GHz processing speed and 8GB RAM. The computer takes less than 1 hour to run the full algorithm.
All simulation results are based on forecasted load as shown Fig.2, with realistic price data taken from [22] and the penalty for each greenhouse gas emission considered from [34]. Fig.3 shows the impact of ω 1 and ω 2 on active power generation of each DG unit at a market price of 24.95 $/ MW. Active power generation from low emission coefficient generators is higher than that from the remaining generators when ω 2 =0.75. This means that DISCO's decision maker gives more priority to emission reduction, and hence generators having high positive impact on emission reduction increase their generation. As shown in Fig.3, active power generation from all Type 2 generators (DG 6-10) increases, as the weight corresponding to emission reduction ω 2 is increased from 0.25 to 0.75. Similarly, as ω 1 increases from 0.25 to 0.75, generators having positive impact on loss reduction increase generation. As can be seen from Fig.3 generation of DG11 is highest when ω 1 =0.75, and decreases with ω 1 decreasing from 0.75 to 0.25, indicating that DG11 has high impact on loss reduction. As the generation of each DG unit is based on its contribution to loss and emission reduction, DISCO's decision maker prioritizes among loss reduction, emission reduction and DISCO's extra benefit.

Impact of ω 1 and ω 2 on a DG unit reactive power price
The impact of ω 1 and ω 2 on Reactive Power Price (RPP) of each DG unit at a market price of 24.95 $/MW is shown in Fig.5. The RPP of each DG unit is based on the contribution of generated reactive power in loss reduction. DG11 has higher impact on loss reduction than the others, and receives more incentives and obtains a higher price for reactive power generation. As DISCO's decision maker increases priority for loss reduction by increasing ω 1 , incentives to each DG unit increase based on contribution to loss reduction leading to an increase in the reactive power price. Fig. 6 shows the variation of losses with respect to ω 1 and ω 2 at a market price of 24.95 $/MW. Increasing ω 1 represents DISCO's decision maker increasing the priority for loss reduction. Hence, losses of network decrease with higher ω 1 value, because of the increase in generation of DG units which have high positive impact on loss reduction.  Fig. 8 shows the variation in DISCO's extra benefit with respect to iterations in the proposed algorithm for different market prices at ω 1 =0.5 and ω 2 =0.5. Zero extra benefit is required to maintain fair competition in a deregulated environment. From Fig. 8, it can be seen that the proposed method provides zero extra benefit at all market prices.

3.7
Comparison in terms of active power loss at ω 1 =0.5 and ω 2 =0.5 The proposed method is compared with the nucleolus theory-based iterative method in [17] in terms of active power loss of the network. As shown in Table 6, at the hours of the day when the market price ( t ) is less than the 'b' coefficients of all DG units the active power losses of the network are the same with both methods. However, at remaining hours of the day active power losses are less in the proposed method than the nucleolus theory-based iterative method.

Comparison in terms of emission
at ω 1 =0.5 and ω 2 =0.5 Table 7 compares emission released from the network with the proposed method and the nucleolus theorybased iterative method. At some hours of the day both methods provide the same emission, because of the wholesale market price ( t ) being less than the 'b' coefficient of all DG units at those hours. As a DG unit will only generate power if the LMP/wholesale market price is greater than its bi coefficient. all DG units are off under such condition and the network is the same as the base case. However, at other hours of the day, the proposed method operates the network with lower emission than the nucleolus-based iterative method.

3.9
Comparison in terms of DG benefit at ω 1 =0.5 and ω 2 =0.5 Fig. 9 compares the DG benefits calculated as the difference between the total amount paid by DISCO to generate power and generation cost, with the proposed and nucleolus theory based iterative methods at a market price of 26.9 $/MW. As seen, the proposed method provides more benefits to generators than the nucleolus theory-based iterative method.

Conclusions
This paper proposes an SVCGT based iterative method for estimating LMPs in DG buses. The LMP value for each DG bus is calculated on the basis of DG's contribution to loss and emissions reduction. DISCO's additional benefit is computed effectively in each iteration using LMP, generation and emission penalty. Two ANN  layers are used to predict load in the next 24 hours based on the load of the previous four hours and at the same hours in the previous two days. The proposed method is developed using the Shapley Value (SV) method, which is monotonous and different to the nucleolus solution concept, and is one of the solutions to the problem of cooperative game theory.
The decision-maker of DISCO can use the proposed method as a control tool for privately owned DG units to help DISCO to operate the network optimally in terms of loss and emission. The proposed method can also help to estimate the status of the network in terms of generation of DG units with controllable additional benefits.
This method can be extended by considering technical objectives such as improvement of reliability and quality of service. The proposed work can also be extended by using the Aumann-Shapley method instead of the Shapley method as the latter requires more computation time. As the integration of DG units into the distribution network is expected to increase in the future, the proposed work can help to resolve issues related to the planning and operation of the distribution network. x : Amount of NO x released from DG unit 'i' (kg/MW); NO Sub x : Amount of No Sub x released based on substation bus load (kg/ MW); P CO 2 : Penalty for the emission of CO 2 in ($/kg); P CO : Penalty for emission of CO in ($/kg); P NO x : Penalty for the emission of NO x in ($/kg); P SO 2 : Penalty for the emission of SO 2 in ($/kg); Ploss t 0 : Base case active power loss at hour t (MW); Ploss t DG : Active power loss of the system with a DG unit at hour t (MW); Ploss t j : System loss with DG units at hour 'h' and iteration 'j'; SO DG i 2 : Amount of SO 2 released from DG unit 'i' (kg/MW); SO Sub 2 : Amount of SO 2 released based on substation bus load (kg/MW); v k (C): Loss(k = l)/emission(k = e) reduction due to coalition C; v k (C − i): Loss(k = l)/emission(k = e) reduction due to coalition C with out DG 'i'; benefit 0 t : Base case DISCO's benefit at hour 't' in ($); benefit j t : DISCO's benefit with DG at hour 't' and iteration 'j' ($); EC 0 t : Green house Gas emission penalty at hour h in the base case in ($); EC DG t : Green house Gas emission penalty at hour h with the DG units in ($); EC j t : Emission penalty for network at hour 't' and iteration 'j' in ($); Emn j t : Emission from network at hour 't' and iteration 'j' in kg; MAPE: Mean Absolute Percentage Error; RMSE: Root Mean Square Error; xEmn(i): Share of DG 'i' in emission reduction (kg); xLoss(i): Share of DG 'i' in loss reduction (MW).