Shapley Value-Based Techno-Economic Framework for Harmonic and Loss Mitigation

This study presents a multi-objective policy for harmonic distortion and loss mitigation in micro-grid and active distribution grid. The implementation of this policy in networks including Distributed Generation (DG), is put forward based on the determination of the locational marginal price of each DG bus, considering their impacts on the loss and harmonic mitigation. In this process, each DG receives an incentive in the form of price increment, based on the Aumann-Shapley value game-theoretic method as an iterative algorithm, considering the effectiveness of reducing loss or harmonic distortion. Furthermore, as a decision-making tool, the operator can change the contribution of the incentives used to mitigate the loss or harmonic distortion and make an appropriate estimation for the next step conditions. The simulations are organized based on real data on the modified Taiwan Power Company grid (TPC). The modified TPC consists of renewable energy sources (RES) and fossil-fuel based DG units. Also, to achieve a realistic framework, generation of RES, loads, and market prices are considered as uncertain parameters. The results demonstrate the competence of the proposed method for the TPC network in terms of triple pricing methods comparison, total loss reduction, harmonic mitigation, and merchandising surplus controlling.


INDEX TERMS
Nowadays one of the most important priorities in the research on distribution networks is loss reduction due to the large number of feeders, connections, and consumers. Distribution network loss is about 10-13% of the total loss. Minimization of loss in distribution networks is studied by many types of research with different solutions [1]. Another critical issue in distribution networks that leads to many problems is harmonic pollution. This issue may result in overheating and tripping on protection devices in transformers, increasing resistance, loss in feeders and conductors, and errors in control and measurement equipment.
The conventional solutions to solve each of these two issues often include costly methods such as adding equipment for harmonic mitigation, e.g. installing Distribution Static Synchronous Compensators (D-STATCOM) [2] or shunt active power filters [3] on distribution networks, and changing the networks configuration for loss reduction [4], [5].
Despite the fact that the presence of PV on the network is sometimes considered as one of the harmonic distortion resources [6], in this paper, we try to utilize PV as a solution to reduce THD level. On the other hand, the main role of power electronic devices, which are utilized as interfaces for RES such as PVs, is providing the power with the required level of power quality. However, in addition to the main role, according to their technical characteristics and the specific decisions of the operators, they can play the auxiliary roles, e.g., harmonic reduction using the virtual harmonic damping impedance concept [7], and selective harmonic elimination pulse-width modulation technique [2], [8].
Using DGs as a problem-solving tool in active distribution networks is a common option in recent research [9], [10]. One of the effective methods that use the DGs as problem-solving tools is pricing method. In this method, there is no need to impose an additional charge, add/remove equipment, or change the network configuration.
The conventional pricing methods in power markets are: uniform [11], zonal [12], marginal loss [13] and locational marginal pricing [14] (LMP). Uniform method considers a single price for all nodes on networks. Therefore, this method cannot be useful for purposeful pricing approaches. Several prices with one price for all nodes in any given zone is called the zonal price [12]. In zonal approach, the market is divided by the congested interconnector. There either needs to be an organized market with a separate price on each side of the interconnector, or there need to be two closely co-operating power exchanges. With marginal loss pricing, transmission losses are priced according to marginal loss factors. The marginal loss factor at a bus is the percentage increase in system losses caused by a small increase in power injection or withdrawal at the bus. Marginal loss factors are always twice average loss factors. However, this pricing method results in an over-collection of loss revenues. Therefore, the fairness of this method in settlement of incentives and benefits may not guaranteed [13]. Locational marginal pricing is a way for wholesale electrical energy prices to reflect the value of electric energy at different locations, accounting for the patterns of load, generation, and the physical limits of the transmission system [15], and the revenue settlement can be performed entirely.
Locational marginal pricing is one of the most effective methods in energy pricing to show the economic impact of each bus on the network. Many of the world's most reliable electricity markets, such as PJM, MISO, and NYISO have employed this pricing method in their transactions [16].
In addition to the common aforementioned problems of harmonic pollution, harmonic flowing in transmission lines and feeders, not only occupies lines capacities and can lead theirs to unexpected congestions, but also the transmission loss is increase. Therefore, the LMP process and the impacts of harmonic pollution are interconnected directly. One of the main contributions of this paper is to reduce the negative effects of harmonic pollution using the purposeful LMP DG Pricing.
In comparison to the uniform [11], [12], zonal [12], and marginal loss [13] methods, LMP can contemplate the technical nature of the network in the economic considerations by utilizing a change in the fair price of each region depending on technical condition. This technical condition can be implemented by the possibility of changing the flow of lines and the possibility of congestion. A comprehensive comparison between these pricing methods are presented in this paper subsequently.
A few number of studies have been conducted in harmonic pricing in power networks. The number of articles is much lower than other technical indices such as loss pricing. In [17], the harmonic pricing is implemented based on the harmonic voltage magnitude. This cost is expressed in the form of a cost function, in which the harmonic voltages cause a negative impact. The first step is harmonic pricing. With harmonic pricing in [18], the customers try to reduce their harmonics. This can be done by substituting nonlinear loads to linear ones or by installing filters at the point of common coupling (PCC). If customers do not perform any of these scenarios, they are forced to pay the penalties. In these papers, the effect of harmonic losses on the pricing procedure has not been clearly and directly seen. Also, the methodologies of these articles are mandatory and did not address the issue of participation in a voluntary manner. In the proposed approach, the ownership of the DGs is considered private; as a result, the process of awarding incentives in the form of price increase for DGs is implemented using a fair game-theoretic method.
In this paper, considering the advantages and shortcomings of the articles mentioned in the literature review, by purposeful DG utilizing in the form of LMP pricing, the total loss and the THD of the network is mitigated simultaneously based on a cooperative game-theoretic approach, without any changes in the network configuration, adding/removing equipment, or budget extension.
The main contributions of this paper are: • Proposing a novel method for harmonic mitigation based on DG pricing, without any changes in the network configuration, adding any new equipment, or budget extension. This procedure for harmonic mitigation is structured by optimal utilization of the power electronic devices connected with DGs as active filters, particularly with PVs. • Comparisons between conventional pricing methods in the facing of loss and harmonic mitigation, and decision making behavior by MS controlling. • Considering the influence of harmonic loss in total loss reduction process in LMP procedure.
The remainder of the paper is structured as follows. Next, the proposed pricing method is presented as four subsections called formulation of harmonic pollution and loss indices, allocation of objectives to DGs, pricing algorithm calculations and backward/forward harmonic power flow. Then, component modeling is outlined in section 3. The uncertainty to the pricing process is investigated in section 4. The simulation results of pricing procedure are addressed in section 5, and the concluding remarks and future recommendations are presented in section 6. A schematic diagram of the pricing procedure, simulation, and the connection between sections of the manuscript is presented in Fig. 1.

II. EXPLANATION OF PROPOSED PRICING STRATEGY
The process of the purposeful locational marginal pricing is presented in this section. In the first part of this section, the formulation of harmonic distortion index and loss as well as the constraints in problem-solving are described. Then, the process of allocating the impact of each DG on mitigating THD and reduction of loss is presented based on the cooperative game-theoretic method. Finally, by combining the relationships between these two sub-sections, the calculations of the iterative algorithm for the locational marginal pricing in distribution networks are presented.

A. FORMULATION OF HARMONIC DISTORTION AND LOSS INDICES
The pricing process based on LMP is studied in this section. The process of LMP pricing based on reducing loss and emission is mentioned in [9] and [19]. However, in this paper for the first time, THD pricing is presented by the LMP method. To maximize the applicability of the proposed framework in the pricing process, in addition to mitigating THD as the first target, loss reduction is considered as the second one.
Considering the objective of THD mitigation and loss reduction in a common framework, while these two objectives are different from each other in terms of kind and quantity, it is necessary to select the similar quantitative indices for these two objectives. This makes it possible to compare and control each one simultaneously or individually. It also makes it possible to simplify the payment process in the next step of proposed algorithm and also allocate financial incentives to DGs. Since the worth of these two objectives can be described by cost, the general objective function is defined by equation (1).
where ''CoL'' represents the net charge for loss and ''CoH'', is the cost that the network must pay due to harmonic negative effects on the equipment, and causing an extra loss in feeders and transmission lines. Harmonic distortions cause remarkable unexpected cost in supply networks as well as at the costumer levels. Estimating the harmonic costs is stand on costs related to harmonic energy losses caused by harmonic pollution.
In addition, the harmonic costs can be defined as the damage of uncontrolled harmonics which is imposed to the system loads and equipment. This extra cost consists of uncontrolled harmonic currents cost, and total harmonic distortion cost caused by uncontrolled harmonic voltage. VOLUME 7, 2019 The effects of these uncontrolled issues are premature aging and derating of equipment. This extra cost is addressed in the form of C(t) in Eq. (2).
Harmonic loss is modeled as an additional resistance which is series with actual resistance. Losses caused by total actual and additional resistances in each line for non-harmonic frequency are equal to losses caused by actual resistance in harmonic pollution condition. Therefore, when optimal power flow (OPF) is performed in fundamental frequency, harmonic effects are not considered in the calculation. To consider the impact of harmonic in OPF, these losses can be modeled in the form of Eq. (3) and (4).
where R both n is the equivalent resistance while considering both power frequency and other harmonic frequencies.
So, for harmonic elimination (or mitigation), it can be possible to spend equal or less than this amount (CoH). The extra harmonic cost proposed by [20], is presented in (5) as follows: Accordingly, for any harmonic distortions, the extra cost called expected value of economic losses due to harmonics, could be calculated according to (5), and the harmonic distortion cost in each duration is determined in the form of Fig. 2 [20].
On the other hand, to calculate the second objective, the cost of loss at the power frequency, according to the FIGURE 2. The harmonic distortion cost Vs duration [20]. energy price in the network (π c ), is obtained from Eq. (6).
Assuming the active and reactive power is extracted from (7): limitations and constraints of this problem are listed in (8) to (12).
• Constraints related to DG power limitations are: • Constraints related to the power factor of the DGs are: • Constraints related to line flow congestion are: • Constraints related to bus voltages are: • Constraints caused by the maximum flow of feeders are:

B. ALLOCATION OF OBJECTIVES TO DGS
After restructuring in power networks, ownership of producers, in particular, the types of DGs, have been privatized [21]. We have tried to establish a fair competitive environment. Existence of fair conditions of the competition leads to the maximum absorption of DG capacity which can be used in decision making process of operators for reducing loss and THD.
The proposed approach is based on Distribution Company (DISCO) performances which are in proportion to the economic signals caused by the loss reduction or THD mitigation. These signals are considered as economic incentives or penalties for DGs. These signals do not interfere with the economic intra-organizational decisions of the DGs, and the economic environment to gain more profits will be quite competitive.
The suggested approach can be described as follows: First, the network is presented without the presence of DGs, and without their effects subsequently. Thus, all the consumers' demand is supplied from the substation. These scenarios are considered as initial and basic conditions of the network. Then the conditions are changed and productions of DGs are injected to the network, and satisfy a part of the system's demand with the substation price. Then, according to the changes in the loss and THD, incentives or penalties are applied to DGs as a pricing procedure.
In the next step, the contribution of each DG to incentives/penalties should be determined. One of the most equitable methods for allocating incentives or penalties to DGs is the application of cooperative game-theoretic methods such as Nucleolus method [22] and Shapley value [23], [24]. In this paper, each DG is considered as a player in the game. The goal of this game is to get more profit for each player participating in supplying the network, as well as reducing loss and mitigating THD. The Shapley value method is employed in this paper. The conventional Shapley value (CSV) has a lot of computational burdens. Therefore, the Aumann method has been used to lighten the high computing burden of CSV method [25]. In the CSV method, for n players, a total set of 2 n − 1 coalitions should be evaluated. Since in this paper, in addition to the allocation of loss, the allocation of harmonic is also considered, the computational burden will become double in the CSV method. Therefore, instead of the CSV method, Aumann-Shapley value method for allocating loss reduction and THD mitigation to each DG is utilized. Equations (13) to (21) explain this approach. N = {1, 2, ..., n} is considered as all player set coalition for n player. In this paper, each DG is considered as a player; S is considered as a subset of N and |s| is the number of players in the S coalition. In this case, Shapley value is in proportion to the loss of player i in the all player set coalition in the equation (13) to (21), according to [26]: If r is the electrical resistance, the real and imaginary parts of branch currents can be written as follows: In this case, the active power loss will be equal to Equation (15).
Subsequently, the unit's Aumann-Shapley partnership index will be as.
Thus, loss reduction can be expressed as (16).
If the allocation is limited to loss, using Aumann method, the active loss allocation to player i from the currents of branches are as Eq. (18).
Nevertheless, in this paper, the main objective is to allocate loss and THD to DGs simultaneously. In this regard, the W(|S|) should be defined precisely.
represents the contribution of each DG to the reduction of loss and the [ν h (s)−ν h (s−d)] represents this contribution to THD mitigation compared to the base status [26]: For THD index, a process is similar to Eq. (17). By analogy, same process is considered in (20): Thus, the contribution of each DG from incentives corresponding to loss reduction ALR d and harmonic mitigation AHM d are determined as follows:

C. PRICING ALGORITHM CALCULATIONS
Considering the possibility of independent decision making by DGs to decrease or increase their generation due to their private ownership, they can affect and change the real-time power flow equations. Therefore, decisions need to be made to predict these conditions. Then the future state of the system is estimated. Under these conditions, DISCO, according to the economic signals derived from the DGs decisions, based on their effects on reducing loss and mitigating THD, executes the following algorithm: Step 1: In this step, because the role of DGs has not been implemented in the network yet, the price of all buses is considered equal to the substation price of the network.
Step 2: According to DG cost function, the generation of DGs are calculated as follows: VOLUME 7, 2019 Step 3: The loss and THD indices are calculated according to the information on the preceding steps and the solution of the power flow equations.
Step 4: The rate of price change, called the incentive, is calculated considering the contribution of each DG in reducing loss and mitigating THD according to Eq. (23).
Step 5: At this step, the LMP of each bus, consisting of the price in the previous iteration (in the first repetition, the price is equivalent to the substation price of the network), and economic incentives, are calculated using (24).
The first terminating constraint (25) is activated when DG generations do not affect loss reduction and THD mitigation. Following the ineffectiveness of increasing DG production in mitigating THD and reducing loss, no incentives are paid and their generations remain constant. If there is no adequate incentive in spite of increasing DG generation, the second terminating constraint will be activated.
In the pricing process, the DISCO's revenue increases by decreasing loss and mitigating THD. This change in the DISCO revenue, compared to the base state is modeled in equation (26) to (28).
The effect of the presence of DG on the changes in the revenue of DISCO is determined by (29) and (30).
In this regard, λ( loss L ) − loss H ) + RTC HM are: The amount added to the DISCO's revenue representing loss reduction and THD mitigation in response to the release of the capacity of the lines respectively. The procedure of the proposed algorithm is depicted in Fig. 3  (RE(i)) are the active and reactive powers at the receiving node of the i th branch as follows: In the above equations, BE is the overall number of buses located beyond the i th branch. Using these equations, the active and reactive power losses in i th branch at h th harmonic are calculated. Therefore, active and reactive powers at the sending node of the i th branch and for h th harmonic can be calculated by the following equations: where r (h) i and x (h) i are the resistance and reactance of the i th branch for h th harmonic.
The current of the i th branch for h th harmonic is given by the following equation:  Therefore, the voltage of the bus at the end of the i th branch is written as follows: To solve the harmonic power flow problem of the pricing procedure, the method proposed in [2] is utilized. Identification of any nodes beyond each branch of the network is realized by a matrix called NBB (Nodes Beyond Branches) and NBN (Nodes Beyond Nodes [2]. The proposed algorithm tries to determine the total generating THD (THD gen ) and the total output of the current (I output ) in each iteration by using a new developed forward/backward harmonic power flow method [2]. The limiting harmonic currents process applied by each DG are calculated according to equations (36) and (37). These values are currents that DGs generate purposefully in order to mitigate the existing harmonics of the network. Then this numerical data is sent to the main pricing algorithm. The flowchart of this procedure is shown in Fig. 3 (highlighted in grey).

III. COMPONENT MODELING A. CONVERTER MODEL
The new wind turbines are usually utilize the AC/DC and DC/AC converters as the interfaces between the generators and the distribution networks [27] and [28]. In the case of higher power requirements, an appropriate option is to use the multi-layer converter topologies [29] and [30]. Similar to the wind turbine, the multi-layer converter topologies with two DC/DC and DC/AC converters [31] share features which make them ideal for the large scale PV system [32].
One of the best harmonic solution options by corrective equipment is active harmonic filters which have excellent cancellation for 2 nd through 50 th harmonic currents, and have no overload problem and, can take advantage of the diversity of loads [33].
As a result, in addition to its power transfer role, the aforementioned PV and wind converters can play the auxiliary role as an active filter for purposes like harmonic mitigation [2]. In this paper, multilayer converters are used for harmonic mitigation.
It is clear that the harmonic mitigation process by employing DGs is based on input current polluted by harmonics and its behavior to eliminate or at least mitigate this undesired factor and deliver unpolluted current. This unpolluted current has only a fundamental component which is recognized by its converter. Therefore, this process can be modeled with three currents called I input as the input current to the device, I gen as the self-generated current of device to mitigate the undesired harmonic component and I output as the desired current which are injected to the network. This modeling is shown in Fig. 4. The current generation process of the device for mitigating the harmonics can be described based on its voltage source converter (VSC).
Considering simultaneous control over pulse width modulation (PWM) and its input DC voltage, it is possible to perform various harmonic elimination schemes like selective harmonic elimination (SHE) which is applied to achieve proper harmonic level by employing precise switching control and algorithm [34]. However, this procedure needs many layers of equipment and this reality leads us to select the multi-layer scheme for SHE process [35]. Therefore, it can be assumed that reducing THD level and SHE can be realized at any desired level corresponds to the cost of the converter [2].

B. NONLINEAR LOAD MODEL
In an industrial distribution system, the main equipment which needs to be dealt with in the harmonic analysis is distribution cables, transformers, nonlinear loads, capacitors and inductors. Instead of using very accurate models, some practical and approximated models for industrial applications are employed in some papers [36], [37]. Therefore, for modeling nonlinear impacts of three-phase rectifier load on the load flow equations, the simulated rectifier is considered as a nonlinear load which is shown in Fig. 5. Thus, for calculating the fundamental or harmonic frequency of bus voltage and line current consequently, the rectifier is modeled as a load. Based on the fundamental and harmonics voltage, the current of the equivalent load is calculated using Eq. (13). The harmonic components can be generated in one phase of VSC, as it is stated in Eq. (14), if the actual DC current is I d [38].

IV. APPLYING UNCERTAINTY TO THE PRICING PROCESS A. DEMAND UNCERTAINTY
The uncertainty modeling methods utilized in electrical energy studies can be classified into two main groups, probabilistic methods such as Monte Carlo [39], [40] and point estimate [41], and non-probabilistic such as robust optimization [42], interval-based analysis [43] and information gap decision theory (IGDT) [43]. Due to the simplicity and the better presentation of random phenomena, and having a closer behavior to the nature of the existing uncertainties in pricing problem, point estimate methods (PEM) is selected for uncertainty modeling in this paper.
This method has split the uncertainty problem into several certain sub-problems. In this way, the nonlinear function of the random problem is converted into k smaller and definite problems only by utilizing k certain variables. In the simplest possible way for this method, two variables can be considered. The first one is higher, and the other one is lower than average. A variable is higher than the mean and another one is lower than that. By this process, for each random variable, the load flow problem is executed twice as a certain procedure.
If it is assumed that X = {x 1 , x 2 , ..., x l , ...x m } is a set of random variables with the mean µ x 1 , standard deviation σ x 1 and stochastic output function of problems is Z = f (x), then for each random variable (x 1 ), the two points which are concentration and weight, corresponding to these variables, are found using the following relationships.
First, the standard locations are found according to (40).
Then, the estimated locations are obtained according to the (41) and (42).
Which the third central moment, based on the probability of occurrence of each state of a random variable, is obtained according to (43).
The stochastic algorithm contained in the flowchart of Fig. 3, which is performed by the PEM method, can be expressed in the following steps: Step I. Determining the number of uncertain input variables Step II. The mean and standard deviation of output variables are considered to be zero.
Step III. The value of the skewness of the variables of step I. is calculated based on the probability of occurrence of each random variable according to (42) and (43).
Step IV. The standard location of random variables is determined according to (40).
Step V. The value of each random variable is determined according to (41).
Step VI. The weighting factor of the random variables is determined according to (42).
Step VII. For each statistical variable, the problem is solved as crisp problems.
Step VIII. Repeat the steps for all uncertain variables of the problem.
The flowchart of the proposed uncertainty algorithm as a stochastic frame is presented in Fig. 3.

B. THE UNCERTAINTY DUE TO WIND POWER
To modeling the uncertainty due to wind power generation, the Weibull probability density function (PDF) is employed to present the stochastic characteristic of wind speed and generation profile, then its PDF and CDF can be formulated as (44) and (45). The estimated Weibull parameters α and β VOLUME 7, 2019 are equal to 12.05 and 2.93 and, the standard error of these parameters are considered equal to 0.39 and 0.21 respectively [20].
Also, the wind power output considering wind velocity is presented as a linear model [44], which is presented as follows: Based on Eq. (46), when the wind speed located between cut-in and rated wind speed, the CDF of wind power output can be expressed as follows: The formulation of the underestimation and overestimation cost of wind power generation is illustrated in [45]. Therefore, these values for each wind generator can be calculated as (48) and (49) respectively [46]: Thus, the total cost of wind power generators can be expressed as follows: where g m is the cost coefficient of the wind generator.

C. THE UNCERTAINTY DUE TO PV
The generation of a PV system stands on environmental parameters: e.g. solar irradiance. The model of this uncertain parameter using a beta PDF is presented in (51) [47].
Based on this distribution function, the output power of the PV system is expressed as a function of irradiation-power curve [48]: In addition, when the solar production exceeds the scheduled one, in order to shift energy purposefully from one time to another, achieve more economic benefits, and obtain a continues THD mitigation process, the PV system is coupled with the energy storage systems.
The state of charge in the electrical energy storage for time interval t in charging/discharging operation mode is a function depending to the charge level of the electrical energy storage at the time interval t − 1 (previous interval), incoming/outgoing power, and the self-discharge rate of the battery. The charge level of the electrical energy storage at each time interval in charging/discharging operation mode is obtained using (53) and (54), respectively.
The total capacity of the battery bank is considered equal to the nominal capacity of all batteries in the system.

V. RESULTS AND DISCUSSIONS
In this section, first, the information of the case study network, 24-hour demands and harmonic loads, and economic characteristics of DGs are presented. Then, the LMP calculation based on DISCO's decisions and the impact of suggested prices on LMP and DG unit generation change is illustrated in Fig. 7 and Fig. 8. In addition, a comparison between pricing methods in the face of loss reduction, THD mitigation, MS control, and line loss reduction is given in Fig. 9 and Fig. 10 respectively.

A. CASE STUDY AND DGS CHARACTERISTICS
In order to evaluate the proposed method, two crisp and stochastic analyses are performed on the 11.4 kV modified Taiwan power company (TPC) 84 bus test network (Fig. 6). The network consists of 11 feeders connected through a bus to the substation. Considering linear and nonlinear loads without the presence of DG, the overall loss and HCOST of this network are equal to 592 kW and 4580 $. Furthermore, this distribution network has twenty-two DG units that can be categorized into five categories. The photovoltaics (PVs), wind turbines, combined cycle gas turbines (CCGT), gas internal combustion engines (G-ICE), diesel internal combustion engines (D-ICE) are presented in table 1. Also, the fuel costs of the two last categories are considered zero because of their renewable nature; the only cost of under/over estimations are applied to the problem. In this paper, the power factor of each unit is considered 0.9 lagging. The data on 24-hour load profile, hourly price, the number of Harmonic Distortions with THD>5%, and the number of interruptions due to harmonic distortions are reported in table 2.

B. LMP CALCULATION BASED ON DISCO'S DECISIONS
The iterative algorithm proposing to solve the LMP pricing problem is a cooperative game in which each DG participant acts as an independent player with the goal of reducing loss and mitigating THD. In this method, the role of the decision maker is presented in the form of parameters w 1 , w 2 . The w 1 represents the importance of reducing loss from the operator's point of view, and similarly, w 2 represents the importance that the operator considers to mitigate the THD of the system. Thus, selecting numbers close to zero for these coefficients means that their importance is very low and the closeness to the number one implies very high importance in the pricing process.
In other words, the DISCO pushes the pricing process towards more loss reduction by increasing the w 1 coefficient, and mitigating THD index leads to a further increment in w 2 . It should be noted that the sum of w 1 and w 2 must be limited to one.

C. IMPACT OF SUGGESTED PRICE ON LMP AND DG UNIT GENERATION CHANGE
The impacts of changing the coefficients on the outputs of each DG unit have been evaluated in Fig. 7 at three suggested prices for the reference bus (20 $/kW, 24 $/kW and 28 $/kW). If the DGs are not active, the same price will be considered for all of them. In other words, at a price of less than 20 $/kW, which is related to the state of the absence of DGs, the DGs have no incentive to activate and consequently do not increase their generations. For this reason, at 20 $/kW, the impact on the pricing process has not been reported and is considered zero in the figure. Therefore, the related results are not shown at this price.
DGs that are far away from the substation, with increasing w 1 and decreasing w 2 , have a more signi*cantimpact on their generations than other DGs due to overall system loss reduction. This trend can be seen in DGs locating on the bus No. 11, 8 and 14 in Fig. 7 (a). As the suggested price increases from 24 $/kW to 28 $/kW, the DG unit's generations increases subsequently, but this procedure does not happen linearly. Furthermore, it is not equal for all these DGs. Therefore, the excessive increase in incentive prices cannot have a significant impact on improving system parameters. This problem is considered as one of the terminating constraints in the proposed iterative algorithm. In this way, the process of applying incentives is adjusted in proportion to the acceptable improvement in each objective. Inversely, if the w 2 increases, then subsequently, the w 1 decreases and the pricing procedure priority changes from loss reduction to THD mitigation.
Finally, if an equal priority is taken for each of the coefficients, the pricing process performs with the same approach to loss reduction and THD mitigation.
The effects of coefficient variations on the LMP of each bus are presented in Fig. 8. In this study, according to the similar procedure presented in Fig. 7, three different types of prices have been evaluated in the process.
Based on the fact that the operational cost of renewable DGs is far lower than the other categories; it is expected that at all suggested prices which are more than substation price, they generate the maximum possible production. However, due to uncertainty and the dependence on natural parameters like wind speed and regarding the low cost of operational cost, the production of renewable DGs are usually smaller than maximum capacity as shown in bus No. 19, 21 in Fig. 7.
On the other hand, by comparison of the results in Fig. 7(a) and Fig. 7(b), it can be seen that despite of the increase in the proposed price from 24 $/kW to 28 $/kW, there is no significant change in wind production. In wind units, due to the insignificant operational cost, the maximum possible power is delivered to the network, regardless of the price difference with the substation price.
Unlike wind power generations, the PV power generations are significantly changed with any change in suggested prices due to using storage. This process is shown on bus No. 18, 20 in Fig. 7. In other words, PV can store the power when the price is low, instead of injecting it to the network. Therefore, price changes have a direct impact on the delivery of PV power to the network. However, if a situation is activated when the PV is scheduled for storage in higher price offering, the networks' objectives can be met by purposeful delivering of PV power to the network. In this study, the maximum storage capacity is considered 20% of the maximum output of the PV.
The economic incentives provided by DISCO directly affect the generation of DGs; in Fig. 8, similar results same as Fig. 7 are expected. The increase in prices indicates the growth of unit generation in accordance with the operator's goals. This procedure reduces the loss and mitigates THD of the system due to the reduction in H COST . Figure 9. compares the conventional pricing methods with the proposed method in the case of loss reduction, THD mitigation, and uncontrolled excess profits of DISCO (merchandising surplus). It is shown that in the uniform pricing method, when the price of all busses is considered to be the same as the substation price, the highest loss and H COST occur in the system. In this method, no decision-making is taken to control the merchandising surplus of DISCO, so the maximum amount of DISCO's uncontrolled surplus is also achieved in this way.

D. COMPARES PRICING METHODS IN THE FACE OF LOSS REDUCTION, THD MITIGATION AND MS CONTROL
In response to generation increment equal to 1 MW, if loss reduction and H COST are investigated, or in other words, the marginal method is used, then compared to a uniform pricing method, loss and H COST are reduced slightly. (The DISCO encourages DGs to generate based on its targets. Nonetheless, since the amounts of these encouragements are not remarkable, DG generation change is slight. Therefore, the DISCO's financial incentives are not impressive and effective). Under these circumstances, considering that DISCO only allocates a small portion of the profit, which is obtained from the rising prices, to DGs; their participation in achieving the goals of the operation are not very significant. Comparing the proposed method with uniform and marginal pricing methods, it can be seen that the amounts of loss and H COST reduction are much lower and, the participation of DGs in meeting the operational goals is increased. Moreover, because of the iterative nature of the proposed method, the uncontrolled merchandising surplus is zero in an iterative descending process.

E. COMPARE PRICING METHODS TO REDUCE LINE LOSS
The line losses must be limited due to the cost imposed on the system. In addition to the conventional loss caused by the current at the fundamental frequency, the flow of other frequencies also leads to further losses in the network. This additional loss in distribution systems can increase power losses up to 20% [49]. Therefore, by releasing the occupied line's capacity and mitigating the harmonic flows through the lines and branches, the total loss of the system declines. Fig. 10 compares the uniform and marginal pricing method with the proposed method in the evaluation of line loss reduction. Since loss are considered in terms of kW, any small change in this figure would be a significant loss reduction.
For example, in the line leading to Bus No. 4 in figure 10(a), the difference between line loss in the uniform and marginal and also between the uniform and proposed pricing methods are 1.786 kW and 4.047 kW respectively. These values for harmonic loss are 3.661 kW and 4.797 kW respectively as shown in Fig. 10(b). According to Fig. 10, the priority of the system management is to reduce harmonics, so the contribution toward harmonic loss reduction is significant. Fig. 10(c), confirms that the proposed method for reducing the loss of lines has been acceptable in comparison with other methods. For example, in the line leading to Bus No. 4, the total loss reductions were 46.49% and 14.81%, compared with the uniform and marginal methods respectively.

VI. CONCLUSION
The efficiency and economic viability of the proposed method have been demonstrated to be superior to other pricing methods through comparative analyses. The more control over the DG units can be implemented by the decision makers' strategy by encouraging them to participate in THD mitigation and loss reduction based on their impacts on satisfying afore-mentioned objectives fairly.
To evaluate the proposed method, and based on the simulation of various pricing methods, it can be concluded that the utilization of LMP method is more efficient than uniform and marginal pricing in terms of loss reduction, THD mitigation and, MS controlling. In comparison with uniform and marginal pricing method, the proposed method decreases the total loss by 1035.41 kW and 633.14 kW respectively. By analogy, the same improvement happens for HCOST reduction and its values decline by 1507.1 $ and 997.9 $ respectively, which verifies the proper performance of the proposed method and its superiority over the other methods in losses and harmonics reduction.
The obtained results from Fig.7 to Fig.9 are depicted that all DG units participate in loss reduction process. However, in the harmonic mitigation procedure, renewable recourses play a crucial role and have a great effect. Further, because of storage system function, PV has the best performance on this issue than other generation units.
Furthermore, it can be seen from Fig. 10 that the proposed method has a significant role in mitigating harmonic losses. These results have also been effective in reducing total loss of the network consequently.
Moreover, the MS in the proposed method is reached to zero approximately, whilst other pricing methods cannot control this parameter suitably.
For future work, DGs can be employed by incentive approach of the proposed method to improve the circumstances of the ancillary services such as reactive power control or power quality indices, e.g. voltage sag/swell, and voltage flicker.